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NOTATION 


A list of the symbols used throughout this document and their definitions is 
provided below for convenience. Parameter values are shown in parentheses. 

Roman Symbols 


a . . . speed of sound 
c . . . airfoil chord 

cp . . . specific heat at constant pressure 
c v . . . specific heat at constant volume 
e . . . internal energy 

i... first grid index of numerical solution 
j ... second grid index of numerical solution 

k... third grid index of numerical solution or thermal conductivity 
k r . . . reduced frequency k r = wc/( 2V) 

Z . . . turbulence model damping function 

„ . time step index of numerical solution or rotational speed (revolutions/sec) 
n... outward unit normal vector 
p . .. pressure 

r . . . radius or cylindrical radial coordinate 
t . . . time 
v . . . velocity 

x ... first Cartesian coordinate 
y . . . second Cartesian coordinate 

z . . . cylindrical axial coordinate or third Cartesian coordinate 
A . . . surface area 

4+ . . . turbulence model parameter 


xi 



B . . . number of propeller blades 

Cp ... power coefficient (Cp = P/pn^D^) 

Ci... thrust coefficient (Cf = T/pn^D^) 

Cep • . . turbulence model parameter 
C kleb • • ’ turbulence model parameter 
C wake • * • turbulence model parameter 

D ... dissipation flux vector, turbulent damping parameter, or diameter 

F . . . flux vector in cylindrical coordinate z direction or turbulence model function 

G . . . flux vector in cylindrical coordinate r direction 
H . . . flux vector in cylindrical coordinate 9 direction 
Hi — total enthalpy 

J . .. advance ratio ( J = U/nD) 

K . .. source term flux vector or turbulence model parameter 
L . . . length 

M . . . Mach number 

P . . . power 

Pr . . . Prandtl number 

B r turbulent • • ♦ turbulent Prandtl number 

Q . . . vector of dependent variables 

R . . . gas constant or residual or maximum radius 

S . . . arc length or pertaining to surface area normal 

T . . . temperature or torque 
U . . . freestream velocity 

V . . . volume 


Greek Symbols 

a . . . time-stepping factor 

2 

e ... modified second-order damping coefficient 
. . . modified fourth-order damping coefficient 
p . . . density 
2 

k ... second-order damping coefficient 
. . . fourth-order damping coefficient 


Xll 



7 . • ■ specific heat ratio 

8 . . . spatial second-order central difference operator 
A . . . blockage factor 

\ v ... second coefficient of viscosity (= — 

ft .. . coefficient of viscosity 

Tf . . . radial transformed variable 

<jj . . . oscillation frequency (normally radians/sec) 

£ . . . axial transformed variable 
( . . . circumferential transformed variable 
v . . . damping factor 
T . . . boundary layer dissipation factor 
A . . . increment of change 


Special Symbols 

V... spatial vector gradient operator 

A . . . spatial forward difference operator (A ~~ ^i) 
y . . . spatial backward difference operator {$ = 4>i ~ l) 

Superscripts 

[ — ] . . . averaged variable 
[ ] . . . dimensional variable 

[ ] . . . implicitly smoothed variable 

[ — *] . . . vector variable 
[ ]* . . . intermediate variable 

[ ] n . . . time step index of variable 

Subscripts 

I]e//ect.»e--' effectivefiowviJ “' 

grid point index of variable 

f 1, • „ . . . laminar flow value 

L J laminar 

[ ]max • • ■ maximum value 


XIII 


[ ] mtn • • . minimum value 
[ ]p . . . related to pressure 
[ ]ps • • • pressure (high pressure) surface 
[ • • • suction (low pressure) surface 

[ . . . total quantity 

[ }z • • • derivative or value with respect to z 

[ ]r • • • derivative or value with respect to r 

[ ]q . . . derivative or value with respect to 0 

[ ]gt a bl e • • • related to stability 
[ kurbulent • • • turbulent flow value 
[ ]oo • * * freestream value 
[ } re f . . . reference value 
[ }kl e l • - • Klebanoff intermittency factor 
[ } wa ke * * * turbulent flow wake parameter 
[ ]2 * . . second-order value 
[ ]4 . . . fourth-order value 


xiv 



ABBREVIATIONS 


A list of the abbreviations used throughout this document and their definitions 

is provided below for convenience. 

ADPAC . . . Advanced Ducted Propfan Analysis Codes 
AOACR... Angle of Attack/Coupled Row aerodynamic analysis code 
ASCII . . . American Standard Code for Information Interchange 
CFL . . . Courant-Freidrichs-Levy number (At/At maa . )S f a j>/ e ) 

MAKEADGRID . . . ADPAC-AOACR multiple-block grid construction program 
ROTGRI D ... Ducted propfan full rotor grid construction program 
SDBLIB . . . Scientific DataBase Library (binary file formats) 

SETUP . ■ ■ ADPAC-AOACR Standard Configuration Setup Program 


xv 




1. SUMMARY 


The primary objective of this study was the development of a time-marching 
three-dimensional Euler/Navier-Stokes aerodynamic analysis to predict steady and 
unsteady compressible transonic flows about ducted and unducted propfan propulsion 
systems employing multiple blade rows. The computer codes resulting from this study 
are referred to as ADPAC-AOACR (Advanced Ducted Propfan Analysis Codes- Angle 
of Attack Coupled Row). This document is the Final Report describing the theoretical 
basis and analytical results from the ADPAC-AOACR codes developed under Task 5 
of NASA Contract NAS3-25270, Unsteady Counterrotating Ducted Propfan Analysis. 

The ADPAC-AOACR program is based on a flexible multiple blocked grid dis- 
cretization scheme permitting coupled 2-D/3-D mesh block solutions with application 
to a wide variety of geometries. For convenience, several standard mesh block struc- 
tures are described for turbomachinery applications. Aerodynamic calculations are 
based on a four-stage Runge-Kutta time-marching finite volume solution technique 
with added numerical dissipation. Steady flow predictions are accelerated by a multi- 
grid procedure. Numerical calculations are compared with experimental data for 
several test cases to demonstrate the utility of this approach for predicting the aero- 
dynamics of modern turbomachinery configurations employing multiple blade rows. 
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2. INTRODUCTION 


Fuel efficient aircraft propulsion systems based on ultra high bypass design tech- 
nologies (propfans, ducted propfans) show great promise for reducing life cycle and 
direct operating costs for both commercial and military transport aircraft. The inte- 
gration of this technology in production aircraft requires extensive aerodynamic and 
structural analyses to verify and optimize designs, and to minimize the deficiencies 
associated with the ultra high bypass concept. Fan noise, blade flutter, and composite 
material behavior are examples of areas which require investigation to lend confidence 
to the ultra high bypass concept. In addition to the above, the overall drag due to the 
large diameter cowl, and potential problems associated with engine placement and 
mounting apparatus, may also be of significant importance. An illustration of the 
aerodynamic characteristics associated with the ultra high bypass propulsion concept 
is given in Fig. 2.1. In these cases, the fan design is typically based on rela- 
tively large, close coupled counterrotating airfoils operating in a high airflow velocity 
environment. Such configurations are highly susceptible to aerodynamic losses at- 
tributable to the interaction between the neighboring blade rows and/or mounting 
hardware, and accurate assessment of the magnitude and influence of such losses is 
crucial in establishing an efficient design. 

Advances in individual component aerodynamic performance in aircraft gas tur- 
bine engine propulsion systems over the past decade has resulted, in part, due to the 
availability of improved computational fluid dynamics (CFD) aerodynamic tools for 
design analysis. Detailed, three-dimensional aerodynamic analyses for isolated turbo- 
machinery blade rows have been presented by Weber and Delaney [1], and Chima [2], 
among others. One deficiency with the isolated blade row approach is that the inter- 
actions between adjacent blade rows in multistage turbomachinery are not properly 
accounted for. These interactions often lead to performance degradation through 
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Figure 2.1: Ultra High Bypass Propulsor aerodynamic characteristics 


4 


incidence swings, hysteresis, and stage mismatching. 

Historically, the prediction of three-dimensional flows through multistage tur- 
bomachinery has been based on one of three solution schemes. These schemes are 
briefly illustrated and described in Figure 2.2. The first scheme involves predicting 

the time-resolved unsteady aerodynamics resulting from the interactions occurring 
between relatively rotating blade rows. Examples of this type of calculation are given 
by Rao and Delaney [39], Jorgensen and Chima [38], and Rai [4]. This approach 
requires ether the simulation of multiple blade passages per blade row, or the incor- 
poration of a phase-lagged boundary condition to account for the differences in spatial 
periodicity for blade rows with dissimilar blade counts. Calculations of this type are 
typically computationally expensive, and are presently impractical for machines with 
more than 2-3 blade rows. 

The second solution technique is based on the average-passage equation system 
developed by Adamczyk [5]. In this approach, separate 3-D solution domains are 
defined for each blade row which encompass the overall domain for the entire tur- 
bomachine. The individual solution domains are specific to a particular blade row, 
although all blade row domains share a common axisymmetric flow. In the solution for 
the flow through a specific blade passage, adjacent blade rows are represented by their 
time and space-averaged blockage, body force, and energy source contributions to the 
overall flow. A correlation model is used to represent the time and space-averaged 
flow fluctuations representing the interactions between blade rows. The advantage of 
the average-passage approach is that the temporally and spatially averaged equation 
system reduce the solution to a steady flow environment; and, within the accuracy 
of the correlation model, the solution is representative of the average aerodynamic 
condition experienced by a given blade row under the influence of all other blade 
rows in the machine. The disadvantage of the average-passage approach is that the 
solution complexity and cost grow rapidly as the number of blade passages increases, 
and the validity of the correlation model is as yet unverified. 

The third approach for the prediction of flow through multistage turbomachin- 
ery is based on the mixing plane concept. A mixing plane is an arbitrarily imposed 
boundary inserted between adjacent blade rows across which the flow is “mixed out” 
circumferentially. This circumferential mixing approximates the time-averaged condi- 
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tion at the mixing plane and allows the aerodynamic solution for each blade passage 
to be performed in a steady flow environment. The mixing plane concept was recently 
applied to realistic turbofan engine configurations by Goyal and Dawes [6]. Flow vari- 
ables on either side of the mixing plane are circumferentially averaged and passed to 
the neighboring blade row as a means of smearing out the circumferential nonunifor- 
mities resulting from dissimilar blade counts. The mixing plane concept is a much 
more cost-effective approach computationally because the flow is steady, and the in- 
dividual blade passage domains are limited to a near-blade region. Unfortunately, the 
accuracy of this approach is clearly questionable under some circumstances because 
of the placement of the mixing plane and the loss of spatial information resulting 
from the circumferential averaging operator. For example, consider a flow field with 
a region of concentrated axial vorticity (such as that resulting from a tip leakage vor- 
tex). Application of the circumferential averaging operator would remove the local 
structure of the vortex in favor of a circumferentially smeared flowfield with equiv- 
alant global properties, resulting in a loss of information concerning the structure of 
the leakage vortex. 

This document contains the Final Report for the ADPAC-AOACR (Advanced 
Ducted Propfan Analysis Codes - Angle of Attack/Coupled Row) 3D Euler/Navier- 
Stokes aerodynamic analysis developed by the Allison Gas Turbine Division of the 
General Motors Corporation under Task V of NASA Contract NAS3-25270. The 
primary objective behind this study was the development of a computational method 
to accurately assess the aerodynamic characteristics of high bypass turbofan engine 
designs employing multiple blade rows. The motivation behind this objective lies in 
the performance potential offered by high bypass ratio turbofan engine designs, and 
the observed effects of aerodynamic interaction occurring between relatively rotating 
adjacent blade rows. 

The ADPAC-AOACR program possesses many features which permit multiple 
blade row solutions using either the time-dependent interaction approach or the mix- 
ing plane concept described above. Average-passage simulations for realistic turbofan 
engine configurations were recently reported under Task 4 of this contract, and further 
details on this approach can be found in Reference [8]. 
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The ADPAC-AOACR program is based on an explicit Runge-Kutta time-marching 
algorithm employing a finite volume, multiple blocked grid discretization. The code 
is constructed such that the algorithm may be applied to multiple blocked grid mesh 
systems with common grid interface (non-overlapping) boundaries. (For a descrip- 
tion of multiple-blocked mesh solution schemes, see Section 3.4.) Several convergence 
enhancements are added for the prediction of steady state flows including local time 
stepping, implicit residual smoothing, and multigrid (See Chapter 3 for further de- 
tails). 

A systematic approach to utilize the ADPAC-AOACR code for the analysis of 
ducted and unducted fan propulsion systems was derived by developing a senes of 
standard blocked mesh configurations for various analytical problems of interest. 
These standard configurations are discussed in detail in the companion report de 
scribing the actual code operation (Computer Program Users Manual) [9]. Several of 
the standard configurations are demonstrated and verified in the discussion of results 

given in Chapter 4.0. 

Separate sections are provided in the chapters which follow to describe the the- 
oretical basis of the ADPAC-AOACR code, and summarize the predicted results and 
verification studies performed to validate the accuracy of the analysis. 

It is worthwhile mentioning that the development and application of the codes 
described in this manual were performed on UNIX-based computers. All files are 
stored in machine-independent format. Small files utilize standard ASCII format, 
while larger files, which benefit from some type of binary storage, are written in 
a machine-independent format through the Scientific DataBase Library (SDBLIB) 
routines [24] The SDBLIB format utilizes machine-dependent input/output routines 
which permit machine independence of the binary data file. The SDBLIB routines 
are under development at the NASA Lewis Research Center. 

Most of the plotting and graphical postprocessing of the solutions was performed 
on graphics workstations. Presently, the PLOTSD [25], SURF [26], and FAST [27] 
graphics software packages developed at the NASA Ames Research Center are being 
extensively nsed for this purpose, and plot output has been tailored for this software. 
In addition, due to the increasing popularity of the PostScript page description Ian- 
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guage, and the variety of devices which can display PostS crip t-based output, a number 
of plotting procedures included in the ADPAC package utilize standard PostScript 
routines. 
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3. 3D EULER/NAVIER-STOKES NUMERICAL ALGORITHM 


In this chapter, the mathematical and computational basis for the ADPAC- 
AOACR time-dependent multiple-grid block Euler /Navier-Stokes analysis is described. 
The definitions of the pertinent variables used in this chapter may be found in the 
Nomenclature. 


3.1 Nondimensionalization 


To simplify the implementation of the numerical solution, all variables are nondi- 
mensionalized by reference values as follows. 

*9 


z — 


J ref 


’ T ~ L 


V = 


Pref 


P = 


ref 


P 

Pref 

f 


ref 


II 

Vr 

vr - » 

v ref 

v e = 

c p 

cv 

k = ■ 

R ref 

Cv — D f ’ 

™ ref 

J 

P 

Pref 

“ L ref 

(Jj ~ 

v ref 



v ref 

k 


z ref 


( 3 . 1 ) 


The reference quantities are defined as follows: 

X , is the maximum diameter of the propfan blade 

re j 

p / is the freestream total pressure 
rrej 

Pref * s f rees ^ ream total density 

a re f is determined from the freestream total conditions 

= y/lPref/Pref 

v f is determined from the freestream total acoustic velocity as 

rej 
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fi re j is determined from the other factors as: 

Pref ~ Pref v ref^ref 

is the freestream thermal conductivity 

R re j is the freestream gas constant 
T re f is the freestream total temperature 


3.2 3-D Navier-Stokes Equations 

The numerical solution procedure is based on an integral representation of the 
strong conservation law form of the Navier-Stokes equations expressed in a cylindrical 
coordinate system. The Euler equations may be derived as a subset of the Navier- 
Stokes equations by neglecting viscous dissipation and thermal conductivity terms 
(i.e. - ft and k = 0). 

Integration of the differential form of the Navier-Stokes equations over a rotating 
finite control volume yields the following equation: 

/ j t (Q)JV + L inv (Q) = f KdV + L ms (Q) (3.2) 


where: 


L inv(Q) ~ [ F inv dA z + ^inv^Ar + ( H{ nv - ru)Q)dAg ] (3.3) 

and: 

Lyis(Q) = J dA [Kis dA z + G vis dA r + H via dA^ (3.4) 

The inviscid (convective) and viscous (diffusive) flux contributions are expressed sep- 
arately by the operators L^ nv and L v ^ a , respectively. 
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The vector of dependent variables Q is defined as: 


Q = 


p 

pVz 

pv r 

P v 6 

■ PH J 


(3.5) 


where the velocity components Vz,Vr, and vq are the absolute velocity components in 
the axial, radial, and circumferential directions relative to the cylindrical coordinate 
system, respectively (see e.g. - Fig. 3.1). The total internal energy is defined 


H = 


(7 ~ l)p ' 2 
The individual flux functions are defined as: 


+ \{ v l + v ? + v|) 


as: 


(3.6) 


F = 
1 tnv 


pv z 


pv r 


P v e 

pv^+p 


pv z v r 


pv z vp 

pv z v r 

o • — 

’ ^inv ~ 

0 

pv r +p 

17. 

’ IJ inv — 

pv T VQ 

PVZVQ 


pv r vg 


(P v p + P ) 

. pvzH . 


. pv r H 


- pvpH . 

r 

0 1 

r ° i 

r 

o 1 


F = 
J vis 


T Z Z 
T zr 

r z0 

L qz J 

F=F(Q), 


n . — 
^ vis ~ 


T rz 
Trr 

T r6 

L qr J 

G = <?(<?), 


IT . 

11 VIS ~ 


T 6z 

T 0r 

T ee 


L 90 J 

H = H(Q) 


(3.7) 


(3.8) 



ADPAC-AOACR Coordinate System Reference 


Cartesian Coordinate 
Reference 


l i 



Figure 3.1: ADPAC-AOACR Cylindrical Coordinate System Reference 
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Finally, the cylindrical coordinate system source term is: 


K = 


0 

pvf+p 


0 


T 96 


(3.10) 


It should be noted that in the numerical algorithm, the radius used in the cylindrical 
source term K must be carefully formulated to guarantee numerical conservation for 
the radial momentum equation. That is, for a uniform flow, the radius in the radial 
momentum equation is chosen such that both sides of the radial momentum equation 
are equal. This ensures that small geometric errors do not corrupt the conservative 
nature of the numerical scheme. The total enthalpy, H, is related to the total energy 
by: 


H — e t + - 

P 


The viscous stress terms may be expressed as: 


(dv z \ 

T zz — 2// ( j + V • V , 

f dvr \ ( dv z \ 

Tzr = >l [(“57 j + (i>r jj ’ 

'■»-*[(;&) + 

^ + A v 


ttt — 2 /i 


'dv r 


V V, 


T rd = 2 P 


:;£)*(&)-(?)]• 


T M = 2tl {r-d6 + r 


VF, 


1 dva dv T \ 

1 ■ I + Av v 

. dT 

qz = V Z T ZZ + V T T ZT + VQT z Q -f k — , 


(3.11) 

(3.12) 

(3.13) 

(3.14) 

(3.15) 

(3.16) 

(3.17) 

(3.18) 
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(3.19) 


dT 


L 

q r = VzTrz + VrTrr + v Q T rO A * Q r ' 


dT 


(3.20) 


qg = v z tq z + v r rg r + vgTgg + k ^ , 
where p is the first coefficient of viscosity, A v is the second coefficient of viscosity, 
and: - a ~ 1 (3.21) 


- du z du r 1 dug u r 
v - v = ~dT + -fr + r-W + T 


The remaining viscous stress terms are defined through the identities: 


Trz — T zr 5 

T 6r = T r6 ’ 
T 0z ~ T z6 » 


(3.22) 

(3.23) 

(3.24) 


3.3 2-D Navier-Stokes Equations with Embedded Blade Row Body 

Forces 

In order to permit a simplified analysis of two-dimensional flows in axisymmetric 
passages with turbomachinery blade rows, a reduced (axisymmetric) form of the 3-D 
Navier-Stokes equations given in the previous section is derived for 2-D flows with 
embedded blade rows represented by axisymmetric blockage and body forces. 

Integration of the differential form of the Navier-Stokes equations for an axisym- 
metric flow yields the following equations: 


f A (A Q)dV + L inv (\Q) = J *SdV + J A KdV + L vis {XQ) 

(3.25) 

where: 



Linv^Q) = J dA Vhnv dAz + MinvdAr 

1 

(3.26) 

and: _ 1 


(3.27) 

LvisW) = J dA [ XF vis dAz + MvisdAr] 

1 
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The term A represents a circumferential blockage factor the value of which may range 
from 0.0 to 1.0. This factor represents the relative circumferential blockage imposed 
by the thickness of the embedded blade row and is calculated as 


A/? 

A = 1.0 (3.28) 

if 

The inviscid (convective) and viscous (diffusive) flux contributions have been ex- 
pressed separately by the operators L^ nv and L v ^ 3 , respectively. 

The vector of dependent variables Q is again defined as: 


P 

pv z 


Q = 


pv r 


P v 6 

■ PH - 


(3.29) 


where the velocity components Vz,vr, and vq are the absolute velocity components in 
the axial, radial, and circumferential directions relative to the cylindrical coordinate 
system, respectively (see e.g. - Fig. 3.1). The total internal energy is defined as: 


e t 


P , 1 ^ 2 . 2 , 2n 
( 7 -l) / > + 2 (V2 + t,r+ ^ ) 


(3.30) 


The vector S contains the body forces and energy sources associated with the 
axisymmetric influence of an embedded blade row. 

The individual flux functions and viscous stresses are as defined in the previous 
section. 

The similarity between the 3-D Navier-Stokes equations and the 2-D axisymmet- 
ric equations described here is obvious. As a result, the numerical implementation 
of each scheme (described in the next section) is essentially the same, and in most 
cases, unless a specific requirement for the 2-D scheme must be described, only the 
3-D numerical scheme details are given. 
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3.4 Numerical Formulation 


The discrete numerical solution is developed from the integral governing equa- 
tions derived in the previous two sections by employing a finite volume solution pro- 
cedure. This procedure closely follows the basic scheme described by Jameson [35]. 
In order to appreciate and utilize the features of the ADPAC-AOACR solution sys- 
tem, the concept of a multiple blocked grid system must be fully understood. It is 
expected that the reader possesses at least some understanding of the concepts of 
computational fluid dynamics (CFD), so the use of a numerical grid to discretize a 
flow domain should not be foreign. Many CFD analyses rely on a single structured 
ordering of grid points upon which the numerical solution is performed (the authors 
are aware of a growing number of unstructured grid solution techniques as well, but 
resist the temptation to mention them in this discussion). Multiple blocked grid sys- 
tems are different only in that several structured grid systems are used in harmony 
to generate the numerical solution. The domain of interest is subdivided into one 
or more structured arrays of hexahedral cells. Each array of cells is referred to as a 
“block”, and the overall scheme is referred to as a multiple blocked mesh solver as 
a result of the ability to manage more than one block. This concept is illustrated 
graphically in two dimensions for the flow through a nozzle in Figures 3. 2-3. 4. 

The grid system in Figure 3.2 employs a single structured ordering, resulting 
in a single computational space to contend with. The mesh system in Figure 3.3 is 
comprised of two, separate structured grid blocks, and consequently, the numerical 
solution consists of two unique computational domains. In theory, the nozzle flowpath 
could be subdivided into any number of domains employing structured grid blocks re- 
sulting in an identical number of computational domains to contend with, as shown in 
the 20 block decomposition illustrated in Figure 3.4. The complicating factor in this 
domain decomposition approach is that the numerical solution must provide a means 
for the isolated computational domains to communicate with each other in order to 
satisfy the conservation laws governing the desired aerodynamic solution. Hence, as 
the number of subdomains used to complete the aerodynamic solution grows larger, 
the number of inter-domain communication paths increases in a corresponding man- 
ner. (It should be noted that this domain decomposition/communication overhead 
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ADPAC-AOACR 2-D Single Block Mesh Structure Illustration 


Physical Domain 



Computational Domain 



Figure 3.2: ADPAC-AOACR 2-D Single Block Mesh Structure Illustration 
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ADPAC-AOACR 2-D Two Block Mesh Structure Illustration 


Physical Domain 



Computational Domain 



Figure 3.3: ADPAC-AOACR 2-D Two Block Mesh Structure Illustration 
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ADPAC-AOACR 2-D Multiple Block Mesh Structure Illustration 

Physical Domain 



Computational Domain 



ja 

Block 


Inter-block communication required 
to couple computational domains 


Figure 3.4: ADPAC-AOACR 2-D Multiple Block Mesh Structure Illustration 
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relationship is also a key concept in parallel processing for large scale computations, 
and thus, the ADPAC-AOACR code appears to be a viable candidate for paralleliza- 
tion via the natural domain decomposition division afforded by the multiple-blocked 
grid data structure.) Clearly, it is often not possible to generate a single structured 
grid to encompass the domain of interest without sacrificing grid quality, and there- 
fore, a multiple blocked grid system has significant advantages. 

The ADPAC-AOACR code was developed to utilize the multiple blocked grid 
concept to full extent by permitting an arbitrary number of structured grid blocks with 
user specifiable communication paths between blocks. The inter-block communication 
paths are implemented as a series of boundary conditions on each block which, in some 
cases, communicate flow information from one block to another. The advantages of 
the multiple-block solution concept are exploited in the calculations presented in 
Chapter 4 as a means of treating complicated geometries with multiple blade rows of 
varying blade number, and to exploit computational enhancements such as multigrid. 

The solution for each mesh block in a multiple-blocked grid is computed iden- 
tically, and therefore the numerical approach is described for a single mesh block. 
In any given mesh block, the numerical grid is used to define a set of hexahedral 
cells, the vertices of which are defined by the eight surrounding mesh points. This 
construction is illustrated in Figure 3.5. 

The cell face surface area normal vector components dAz, dAr, and dAy are 
calculated using the cross product of the diagonals defined by the four vertices of 
the given face, and the cell volume is determined by a procedure outlined by Hung 
and Kordulla [33] for generalized nonorthogonal cells. The integral relations ex- 
pressed by the governing equations are determined for each cell by approximating the 
area-integrated convective and diffusive fluxes with a representative value along each 
cell face, and by approximating the volume-integrated terms with a representative 
cell volume weighted value. The discrete numerical approximation to the governing 
equation then becomes 


Vol 


At 


(Finv(Q)i+l,j,k Finv(Q)i,j,k 


(3.31) 


+Ginv(Q)i,j+l,k 


(*inv 
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Grid Point 

Grid Line 

Hidden Grid Line 

Surface Area Normal Vector 

Cell Face Diagonal Vector 



*•* Volume V 

Figure 3.5: Three-dimensional finite volume cell 
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+ H inv(Q\,j,k+\ ~ H inv(Q)i,j,k 

+Fvis(Q)i+l,j,k ~ Fvis(Q\,j,k 

+G v is(Q)i,j + l,k ~ & vis(Q)i,j,k 
+tfvis(0)i,j,fc+ 1 _ H vis(Q)i,j,k) 

+VolK + Di^ k (Q) 

Here, t,j,fc represents the local cell indices in the structured cell array, Vol is 
the local cell volume, At is the calculation time interval, and D { j k is an artificial 
numerical dissipation function which is added to the governing equations to aid nu- 
merical stability, and to eliminate spurious numerical oscillations in the vicinity of 
flow discontinuities such as shock waves. Following the algorithm defined by Jame- 
son [35], it is convenient to store the flow variables as a representative value for the 
interior of each cell, and thus the scheme is referred to as cell-centered. The discrete 
convective fluxes are constructed by using a representative value of the flow variables 
Q which is determined by an algebraic average of the values of Q in the cells lying on 
either side of the local cell face. Viscous stress terms and thermal conduction terms 
are constructed by applying a generalized coordinate transformation to the governing 

equations as follows: 

4 = £(z,r,0), 77 = r?(z,r,0), C — C( z i r i^) ( 3 ’ 32 ) 


The chain rule may then be used to expand the various derivatives in the viscous 


stresses as: 
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(3.33) 

(3.34) 

(3.35) 


The transformed derivatives may now be easily calculated by differencing the variables 
in computational space (1 corresponds to the f direction, j corresponds to the >) di- 
rection, and k corresponds to the ( direction), and utilizing the appropriate identities 
for the metric differences (see e.g. [32]). 
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3.5 Runge-Kutta Time Integration 


The time-stepping scheme used to advance the discrete numerical representation 
of the governing equations is a multistage Runge-Kutta integration. An m stage 
Runge-Kutta integration for the discretized equations is expressed as: 

Q 1 =Q n -a 1 At{L(Q n ) + D(Q n )}, 

Q 2 = Q n - <* 2 At[X(Qi) + D(Q n )}, 

Qz = Q n - a 3At[I(C?2) + D(Q n )]i 

Q A = Q n - a 4 At[X(Q 3 ) + D(Q n )], 

Qm — Q n ~ a m M[L{Q m _i) + D(Q n )}, 

Q n+1 = Qm (3-36) 

where: 

HQ) = Hnv(Q) - < 3 - 37 ) 

For simplicity, viscous flux contributions to the discretized equations are only cal- 
culated for the first stage, and the values are frozen for the remaining stages. This 
reduces the overall computational effort and does not appear to significantly alter the 
solution for unsteady flows (there is no difference for steady state flows). It is also 
generally not necessary to recompute the added numerical dissipation terms during 
each stage. In this study, three different multistage Runge-Kutta schemes (2 four- 
stage schemes, and 1 five-stage scheme) were investigated. The differences between 
the schemes will be described in more detail in a later section. 

A linear stability analysis of the four and five-stage Runge-Kutta time-stepping 
schemes utilized during this study indicate that the schemes are stable for all calcu- 
lation time increments St which satisfy the stability criteria CFL < 2\/2 (four stage) 
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or CFL < 3.75 (five stage). Based on convection constraints alone, the CFL number 
may be defined in a one-dimensional manner as: 


CFL = 


At 




(3.38) 


In practice, the calculation time interval must also include restrictions resulting from 
diffusion phenomena. The time step used in the numerical calculation results from 
both convective and diffusive considerations and is calculated as: 


At = CFL 


H ) 

A; -I- A j + + I/; + Uj + V ) (j y 


(3.39) 


where the convective and diffusive coordinate wave speeds (A and v, repsectively) are 
defined as: 

A • = Vol/{V -Si + a) (3.40) 

pVofl 

Vi ~ C At (S*)p 

The factor C aj is a “safety factor” of sorts, which must be imposed as a result of 
the limitations of the linear stability constraints for a set of equations which are truly 
nonlinear. This factor was determined through numerical experimentation and is 
normally chosen to be 4.0. 

For steady state flow calculations, an acceleration technique known as local time 
stepping is used to enhance convergence to the steady-state solution. Local time step- 
ping utilizes the maximum allowable time increment at each point during the course 
of the solution. While this destroys the physical nature of the transient solution, the 
steady-state solution is unaffected and can be obtained in fewer iterations of the time- 
stepping scheme. For unsteady flow calculations, of course, a uniform value of the 
time step At must be used at every grid point to maintain the time-accuracy of the 
solution. Other convergence enhancements such as implicit residual smoothing and 
multigrid (described in later sections) are also applied for steady flow calculations. 
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3.0 Fluid Properties 


The working fluid is assumed to be air acting as a perfect gas, thus the ideal gas 
equation of state has been used. Fluid properties such as specific heats, specific heat 
ratio, and Prandtl number are assumed to be constant. The fluid viscosity is either 
specified as a constant, or derived from the Sutherland (see e.g. [32]) formula: 


p = c x 


3 

(gg 

t + c 2 


The so-called second coefficient of viscosity A r is fixed according to: 


(3.42) 


A 



(3.43) 


The thermal conductivity is determined from the viscosity and the definition of the 
Prandtl number as: 


k = 


c pP 

Pr 


(3.44) 


3.7 Dissipation Function 

In order to prevent odd-even decoupling of the numerical solution, nonphysical 
oscillations near shock waves, and to obtain rapid convergence for steady state solu- 
tions, artificial dissipative terms are added to the discrete numerical representation 
of the governing equations. The added dissipation model is based on the combined 
works of Jameson et al. [35], Martinelli [10], and Swanson et al. [36]. A blend of 
fourth and second differences is used to provide a third order background dissipation 
in smooth flow regions and first order dissipation near discontinuities. The discrete 
equation dissipative function is given by: 

= ( D i ~ D i + D j ~ D j + D t~ D \ (3-45) 
The second and fourth order dissipation operators are determined by 

D l Q w = A* (3.46) 
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( 3 . 47 ) 


D (Q id* = J,k ] A < A < Qi ’ j ’ k 

wh ete Ac and V £ are forward and backward difference operators in the £ direction 
In order to avoid excessively large levels of dissipation for cells w,th large aspect 
ratios, and to maintain the damping properties of the scheme, a vanable scalmg of 
the dissipative terms is employed which is an extension of the two dimens, ona. scheme 
given by Martinelli [ 10 ]. The scaling factor is defined as a funct.on of e spe 
radius of the Jacobian matrices associated with the £, ,, and ( directions and P~v.de, 
a scaling mechanism for varying cell aspect ratios through the follow, ng scheme. 

( 3 . 48 ) 




The function $ controls the relative importance of dissipation in the three coordinate 
directions as: 


$ v _ = 1 + max | (tt-, ) s ’ 


( 3 . 49 ) 


The directional eigenvalue scaling functions are defined by: 

(*f ) i+ l j,k = + 

(A,) «+^j ,k = u i+y,k {Sv) i+y,>‘ + dS '' ) ' + 2 j ' t 
= v i+y /Wi+y* + c(S() *+2J,*= 

The use of the maximum function in the definition of * is important for grids where 
A„/Ar and Ac /Ac are very large and of the same order of magnitude. In tk. case^ 
if th L ratios are summed rather than taking the maximum, the d.ss.pat.on can 
become too large, resulting in degraded solution accuracy and poor 
Because three-dimensional solution grids tend to exhibit large vanat.ons ,n the 
aspect ratio, there is less freedom in the choice of the parameter a for tbs scheme, 
and a value of 0.5 was found to provide a robust scheme. 


( 3 . 50 ) 

( 3 . 51 ) 

( 3 . 52 ) 
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The coefficients in the dissipation operator use the solution pressure as a sensor 
for the presence of shock waves in the solution and are defined as: 


£ - + l j >k = 

v . _ ~ + Pi+l,j,k)\ 

l '’ > ' ( p i~ 1 ,j,k + 2 p i,j,k + Pi + 1 ) 

e4 |l -l = n»a*(0,K 4 -e? I ) 

9 4 

where k are user-defined constants. Typical values for these variables are 

2 = I 4 = i_ 

2 64 

The dissipation operators in the 77 and £ directions are defined in a similar manner. 


(3.53) 

(3.54) 

(3.55) 


(3.56) 


3.8 Implicit Residual Smoothing 

The stability range of the basic time-stepping scheme can be extended using 
implicit smoothing of the residuals. This technique was described by Hollanders 
et al. [37] for the Lax-WendrofF scheme and later developed by Jameson [35] for 
the Runge-Kutta scheme. Since an unsteady flow calculation for a given geome- 
try and grid is likely to be computationally more expensive than a similar steady 
flow calculation, it would be advantageous to utilize this acceleration technique for 
time-dependent flow calculations as well. In recent calculations for two dimensional 
unsteady flows, Jorgensen and Chima [38] demonstrated that a variant of the im- 
plicit residual smoothing technique could be incorporated into a time-accurate ex- 
plicit method to permit the use of larger calculation time increments without ad- 
versely affecting the results of the unsteady calculation. The implementation of this 
residual smoothing scheme reduced the CPU time for their calculation by a factor 
of five. This so-called time-accurate implicit residual smoothing operator was then 
also demonstrated by Rao and Delaney [39] for a similar two-dimensional unsteady 
calculation. Although this “time-accurate” implicit residual smoothing scheme is not 
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developed theoretically to accurately provide the unsteady solution, it can be demon- 
strated that errors introduced through this residual smoothing process are very local 
in nature, and are generally not greater than the discretization error. 

The standard implicit residual smoothing operator can be written as: 

(1 - A^ y^)(l - e 7 1 A 77 Vt?)(1 ~ S7()Rm = Rm (3.57) 

where the residual Rm is defined as: 

R m = amyiQm - D m ), rn = 1 ,mstages (3.58) 

for each of the m stages in the Runge-Kutta multistage scheme. Here Qm is the sum 
of the convective and diffusive terms, Dm the total dissipation at stage m, and Rm 
the final (smoothed) residual at stage m. 

The smoothing reduction is applied sequentially in each coordinate direction as. 

Rm = (1 - 

Rm — (1 — e V ^ V V17) ^ Rm 

RT = ( 1-«C V C r^m 

Rm = RT ( 3 ‘ 59 ) 

where each of the first three steps above requires the inversion of a scalar tridiag- 
onal matrix. The application of the residual smoothing operator varies with the 
type of Runge-Kutta time marching scheme selected. The full four and five stage 
time-marching schemes utilize residual smoothing at each stage of the Runge-Kutta 
integration. The reduced four stage scheme employs residual smoothing at the second 
and fourth stages only. 

The use of constant coefficients (e) in the implicit treatment has proven to be 
useful, even for meshes with high aspect ratio cells, provided additional support such 
as enthalpy damping (see [35]) is introduced. Unfortunately, the use of enthalpy 
damping, which assumes a constant total enthalpy throughout the flowfield, cannot 
be used for an unsteady flow, and many steady flows where the total enthalpy may 
vary. It has been shown that the need for enthalpy damping can be eliminated by using 
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variable coefficients in the implicit treatment which account for the variation of the 
cell aspect ratio. Martinelli [10] derived a functional form for the variable coefficients 
for two-dimensional flows which are functions of characteristic wave speeds. In this 
study, the three-dimensional extension described by Radespiel et al. [36] is utilized, 
and is expressed as: 



CFL 1 - t - ma *< r ^’« ) ,2 
e t - max | 0, 4 1 CFLmax 1 + mai(r^r^) 1 

1 CFL 

ejj max | , ^ Lniax 1 + m«r(r^r^) 


e> = max 


0,t1 


CFL 1 + maa ( r i^ r ^ 2 


(3.60) 


(3.61) 


(3.62) 


4 CFLmax 1 + max(r rj ^r^) 

CFL represents the local value of the CFL number based on the calculation time in- 
crement At, and CFLmax represents the maximum stable value of the CFL number 
permitted by the unmodified scheme (normally, in practice, this is chosen as 2.5 for 
a four stage scheme and 3.5 for a five stage scheme, although linear stability analy- 
sis suggests that 2^2, and 3.75 are the theoretical limits for the four and five stage 
schemes, respectively). From this formulation it is obvious then that the residual 
smoothing operator is only applied in those regions where the local CFL number 
exceeds the stability-limited value. In this approach, the residual operator coefficient 
becomes zero at points where the local CFL number is less than that required by 


stability, and the influence of the smoothing is only locally applied to those regions 
exceeding the stability limit. Practical experience involving unsteady flow calcula- 
tions suggests that for a constant time increment, the majority of the flowfield utilizes 
CFL numbers less than the stability-limited value to maintain a reasonable level of 
accuracy. Local smoothing is therefore typically required only in regions of small grid 
spacing, where the stability-limited time step is very small. Numerical tests both with 
and without the time-accurate implicit residual smoothing operator for the flows of 
interest in this study were found to produce essentially identical results, while the 
time-accurate residual smoothing resulted in a decrease in CPU time by a factor of 
2-3. In practice, the actual limit on the calculation CFL number were determined to 
be roughly twice the values specified for C F Lmaxt above. 
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3.9 Turbulence Model 


As a result of computer limitations regarding storage and execution speed, the 
effects of turbulence are introduced through an appropriate turbulence model and 
solutions are performed on a numerical grid designed to capture the macroscopic 
(rather than the microscopic) behavior of the flow. A relatively standard version of 
the Baldwin-Lomax [34] turbulence model was adopted for this analysis. This model 
is computationally efficient, and has been successfully applied to a wide range of 
geometries and flow conditions. 

The effects of turbulence are introduced into the numerical scheme by utiliz- 
ing the Boussinesq approximation (see e.g. [32]), resulting in an effective calculation 
viscosity defined as: 


P" effective P laminar P turbulent (3.63) 

The simulation is therefore performed using an effective viscosity which combines the 
effects of the physical (laminar) viscosity and the effects of turbulence through the 
turbulence model and the turbulent viscosity ^turbulent' 

The Baldwin-Lomax model specifies that the turbulent viscosity be based on an 
inner and outer layer of the boundary layer flow region as: 

{ if 1 turbulent) inner » V — Vcrossover 

(3.64) 

\Hurbulent> outer ■> V > Vcrossover 

where y is the normal distance to the nearest wall, and ycrossover is the smallest 
value of y at which values from the inner and outer models are equal. The inner and 
outer model turbulent viscosities are defined as: 

(Vturb)inner = P* 2 M ( 3 - 65 ) 

(A turb) outer = ^^ c PP^wake^kleby (3.66) 

Here, the term l is the Van Driest damping factor 

l = At/( 1 — e ^~ y ! A )) (3.67) 
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w is the vorticity magnitude, F wake is defined as: 


*wake ~ VmaxFmax 

where the quantities ymax , F m ax are determined from the function 

F(y) = yH[l- e (-S/ + M + )] 

The term y+ is defined as 


(o.OOj 


(3.69) 


y 




p\w\ 


Haminar (3 ' 70) 

The quantity F MAX is the maximum value of F(y) that occurs in a profile, and 
y MAX 1S the value °f V at which it occurs. The determination of F MAX and y^AX 
is perhaps the most difficult aspect of this model for three-dimensional flows. The 
profile of F(y) versus y can have several local maximums, and it is often difficult to 
establish which values should be used. In this case, F MAX is taken as the maximum 

value of F(y) between a y+ value of 350.0 and 1000.0. The function F Ueb is the 
Klebanoff intermittency factor given by 

F kUb (y) = [1 + 


Vmax 

and the remainder of the terms are constants defined as 


(3.71) 


A + = 26, 


Cep — 1 . 6 , 

C kleb = °- 3 > 
k 0.4, 


K - 0.0168 (3.72) 

In practice, the turbulent viscosity is limited such that it never exceeds 1000.0 times 
the laminar viscosity. 

The turbulent flow thermal conductivity term is also treated as the combination 
of a laminar and turbulent quantity as: 


^ effective ^laminar + ^turbulent 


(3.73) 
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For turbulent flows, the turbulent thermal conductivity k turbu i ent is determined from 
a turbulent Prandtl number Pr tur bul en t such that 


P'turbuUnt - k turMent 


(3.74) 


The turbulent Prandtl number is normally chosen to have a value of 0.9. 

In order to properly utilise this turbulence model, a fairly large number of gnd 
cells must be present in the boundary layer flow region, and, perhaps o greater im- 
portance, the spacing of the first grid cell off of a wall should be smaU enoug o 
accurately account for the inner “law of the wall” turbulent boundary layer p 
region. Unfortunately, this constraint is typically not satisfied due to fl-l-mduced 
problems or excessive computational costs, especially for time-dependent flow cak- 
lations. A convenient technique to suppress this problem is the use of wall functmn 
to replace the inner turbulent model function, and solve for the flow on a some 
coarser grid. The wall function technique has been successful y 
computer programs based on a numerical algorithm similar to the ADPA f A0 
solution scheme. Unfortunately, this technique has no. been implemented or tested 
for the current application. It appears that the use of wall functions is a reasonable 

area for future research. ... 

Practical applications of the Baldwin-Lomax model for three-dimensional viscous 

flow must be made with the limitations of the model in mind. The win omax 

mode, was designed for the prediction of wall bounded turbulent shear 

not likely to be well suited for flows with massive separations or large vortical struc- 

lef There are, unfortunately, a number of applications for ducted and unducted 

propfans where this model is likely to be invalid. This is also likely to be an area 

requiring improvement in the future. 

3.10 Multigrid Convergence Acceleration 

Multigrid (not to be confused with a multiple blocked grid!) is a numerical so- 
lution technique which attempts to accelerate the convergence of an iterative process 
(such as a steady flow prediction using a time-marching scheme) by computing ^ 
ructions to the solution on coarser meshes and propagating these changes to 
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mesh through interpolation. This operation may be recursively applied to several 
coarsenings of the original mesh to effectively enhance the overall convergence. In 
the present multigrid application, coarse meshes are derived from the preceding finer 
mesh by eliminating every other mesh line in each coordinate direction as shown in 
Figure 3.6. As a result, the number of multigrid levels (coarse mesh divisions) is 
controlled by the mesh size, and, in the case of the ADPAC-AOACR code, also by 
the indices of the embedded mesh boundaries (such as blade leading and trailing 
edges, etc.) (see Figure 3.6). These restrictions suggest that mesh blocks should be 
constructed such that the internal boundaries and overall size coincide with numbers 
which are compatible with the multigrid solution procedure (i.e., the mesh size should 
be 1 greater than any number which can be divided by 2 several times and remain 
whole numbers: e.g. 9, 17, 33, 65 etc.) 

The multigrid procedure is applied in a V-cycle as shown in Figure 3.7, whereby 
the fine mesh solution is initially “injected” into the next coarser mesh, the appro- 
priate forcing functions are then calculated based on the differences between the 
calculated coarse mesh residual and the residual which results from a summation of 
the fine mesh residuals for the coarse mesh cell, and the solution is advanced on the 
coarse mesh. This sequence is repeated on each successively coarser mesh until the 
coarsest mesh is reached. At this point, the correction to the solution ( Q™ "t " } — Q 1 } . ) 
is interpolated to the next finer mesh, a new solution is defined on that mesh, and the 
interpolation of corrections is applied sequentially until the finest mesh is reached. 
Following a concept suggested by Swanson et al. [36], it is sometimes desirable to 
smooth the final corrections on the finest mesh to reduce the effects of oscillations 
induced by the interpolation process. A constant coefficient implementation of the 
implicit residual smoothing scheme described in Section 3.5 is used for this purpose. 
The value of the smoothing constant is normally taken to be 0.2. 

A second multigrid concept which should be discussed is the so-called “full” 
multigrid startup procedure. The “full” multigrid method is used to initialize a solu- 
tion by first computing the flow on a coarse mesh, performing several time-marching 
iterations on that mesh (which, by the way could be multigrid iterations if succes- 
sively coarser meshes are available), and then interpolating the solution at that point 
to the next finer mesh, and repeating the entire process until the finest mesh level 
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Multigrid Algorithm Mesh Level Description 



Grid lines defining mesh boundaries and internal boundaries (blade leading edges, trailing edges, etc.) 
must be consistent with the mesh coarsening process (cannot remove a mesh line defining a boundan 

for the given coordinate direction) 

Figure 3.6: Multigrid Mesh Coarsening Strategy and Mesh Index Relation 


Multigrid V-Cycle Strategy 

Advance Solution on Current Mesh 
Update Solution with Corrections 


Mesh Level 



Figure 3.7: Multigrid V Cycle Strategy 
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is reached. The intent here is to generate a reasonably approximate solution on 
the coarser meshes before undergoing the expense of the fine mesh multigrid cycles. 
Again, the “full” multigrid technique only applies to starting up a solution. 

3.11 Single Block Boundary Conditions 

In this section, the various boundary conditions which are available as part of 
the ADPAC-AOACR code pertaining to a single mesh block are described. Before 
describing the individual boundary conditions, it may be useful to describe how the 
boundary conditions are imposed in the discrete numerical solution. Finite volume so- 
lution algorithms such as the ADPAC-AOACR program typically employ the concept 
of a phantom cell to impose boundary conditions on the external faces of a particular 
mesh block. This concept is illustrated graphically for a 2-D mesh representation in 
Figure 3.8. 

A phantom cell is a fictitious neighboring cell located outside the extent of a 
mesh which is utilized in the application of boundary conditions on the outer bound- 
aries of a mesh block. Since flow variables cannot be directly specified at a mesh 
surface in a finite volume solution (the flow variables are calculated and stored at 
cell centers), the boundary data specified in the phantom cells are utilized to control 
the flux condition at the cell faces of the outer boundary of the mesh block, and, in 
turn, satisfy a particular boundary condition. All ADPAC-AOACR boundary con- 
dition specifications provide data values for phantom cells to implement a particular 
mathematical boundary condition on the mesh. Another advantage of the phantom 
cell approach is that it permits unmodified application of the interior point scheme 
at near boundary cells. 

Inflow and exit boundary conditions are applied numerically using characteristic 
theory. A one-dimensional isentropic system of equations is utilized to derive the 
following characteristic equations at an axial inflow/outflow boundary: 


dc~ , ,dc~ n 

at (vz a) a* = °’ 

(3.75) 

dc+ , ,dc+ n 

at +( ^ + fl) a, =0 

(3.76) 
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2-D Mesh Block Phantom Cell Representation 


• Grid Point 

Grid Line 

Mesh Block Boundary 
— — Phantom Cell Representation 



"Comer" phantom cells cannot be controlled 
through boundary conditions, but must be updated 
to accurately compute grid point averaged values 


Figure 3.8: 2-D Mesh Block Phantom Cell Representation 
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where: 


2 a 


(3.77) 


C -v z - 


7 


- 1 ’ 


C = t>z + 


2a 


For subsonic normal inflow, the upstream running invariant C~ is extrapolated 
to the inlet, and along with the equation of state, specified total pressure, total tem- 
perature, and flow angle (used to specify the angle of attack), the flow variables at the 
boundary may be determined. Inlet boundaries also require the specification of a flow 
direction. This is accomplished in one of two different manners. For turbomachinery 
based flow calculations, it is typical to specify a radial and circumferential (swirl) 
flow angle. For external flow calculations (such as angle of attack flow analysis) a 
Cartesian based flow angle is prescribed. It should be mentioned that the effective 
inflow angle in this case may vary for a given block as it rotates about the axis, and 
therefore the inflow angle is actually a function of circumferential position, 9. 

Outflow boundaries require a specification of the exit static pressure at either the 
bottom or top of the exit plane. The remaining pressures along the outflow boundary 
are calculated by integrating the radial momentum equation: 


dp = ^_ | (3.78) 

dr r 


In this case, the downstream running invariant (7+ is used to update the phantom 
cells at the exit boundary. 

Far- field boundaries also use the characteristic technique based on whether the 
local flow normal to the boundary passes into or out of the domain. 

All solid surfaces (hub, cowl, airfoils) must satisfy either flow tangency for inviscid 

flow: 

V -n = 0 ( 3 - 79 ) 


or no slip for viscous flows: 

Vz = 0, v r = 0, VQ = tuj (3.80) 

In both cases, no convective flux through the boundary (an impermeable surface) is 
permitted. The phantom cell velocity components are thus constructed to ensure that 
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the cell face average velocities used in the convective flux calculation are identically 
zero. The phantom cell pressure is simply extrapolated based on the boundary layer 
flow concept dp/dn = 0 (although for inviscid flows a normal momentum equation is 
likely to be more accurate). In the present applications, for inviscid flows, extrapo- 
lation was found to be the most effective technique based on rapid convergence and 
solution accuracy. The phantom cell density or temperature is imposed by assuming 
either an adiabatic surface dT/dn = 0 or a specified surface temperature, which sug- 
gests that the phantom cell temperature must be properly constructed to satisfy the 
appropriate average temperature along the surface. 

Some final comments concerning boundaries are in order at this point. Artificial 
damping is applied at the block boundaries by prescribing zero dissipation flux along 
block boundaries to maintain the global conservative nature of the solution for each 
mesh block. Fourth order dissipation fluxes at near boundary cells are computed using 
a modified one-sided differencing scheme. Implicit residual smoothing is applied at the 
block boundary by imposing a zero residual gradient (i.e. ( dR/dz ) = 0.0) condition 
at the boundary. 


3.12 Multiple-Block Coupling Boundary Conditions 

For the multiple-block scheme, the solution is performed on a single grid block 
at a time. Special boundary conditions along block boundaries are therefore required 
to provide some transport of information between blocks. This transport is provided 
through one of four types of procedures. Each procedure applies to a different type 
of mesh construction and flow environment, and details of each approach are given 
in the paragraphs below. 

For neighboring mesh blocks which have coincident mesh points along the in- 
terface separating the two blocks, a simple direct specification of the phantom cell 
data based on the near boundary cell data from the neighboring block has been used 
successfully. This concept is illustrated graphically in Figure 3.9. Each phantom 

cell in the block of interest has a direct correspondance with a near boundary cell in 
the neighboring mesh block, and the block coupling is achieved numerically by simply 
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Contiguous Mesh Block Interface 
Boundary Coupling Scheme 



Phantom cell data values for mesh block# 1 
are determined by corresponding near— boundary 
data in mesh block # 2 



Block Boundary 
Mesh Line 

Phantom Cell Representation 
Phantom Cell Data Path 


Figure 3.9: ADPAC-AOACR Contiguous Mesh Block Coupling Scheme 
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^signing the vnlne of the corresponding cell in the neighboring block to the phantom 
ce o t e block of interest. This procedure essentially duplicates the interior point 
solution scheme for the near boundary cells, and uniformly enforces the conservation 
principles implied by the governing equations. 

For mesh block interfaces which do not have coincident points, but which define 
essentially the same surface, it is necessary to perform some type of interpolation 
to determine representative values for the block-specific phantom cells absed on the 
near boundary interface data from the neighboring mesh block. A relatively simple 
interpolation scheme based on an electrical resistance analogy was developed for this 
purpose, and is illustrated graphically for a non-contiguous mesh block interface as 
s own in Figure 3.10. To determine the phantom cell data for a non-contiguous 

mesh block interface, the center of the nearest cell face for the particular phantom cell 
is determined by averaging the coordinates of the four mesh points which determine 
e cell face. A search is then performed in the neighboring mesh block along the 
same interface to determine the cell face whose center lies closest to the phantom cell 
ace center of interest. It is assumed that both mesh blocks share a common surface 
at the interface of interest, in spite of the fact that this surface is defined by different 
points. This assumption thus requires that the block boundaries do not overlap, and 
t at the physical extent of the boundaries are approximately equal. Once the newest 
ce in the neighboring mesh block has been located, the nearest set of four cells 
“ thC nei S hborin 6 mesh is determined, and an electrical circuit analogy is used to 
interpolate the phantom cell boundary data from the four cells of neighboring mesh 
block data. The circuit analogy simply sets a resistance for each path of neighboring 
cell data to the phantom cell center based on the physical distance between the centers 
o each cell face (again, this concept is illustrated in Figure 3.10. The four neighboring 
cell near boundary data are then treated as voltages, and the phantom cell data are 
etermined by calculating the appropriate “voltage” potential at the phantom cell 
ace center based on the neighboring cell face center data and the resistance along 
each path. This scheme has the advantages of being simple to program, is extremely 
stable, and identically duplicates the contiguous mesh block interface phantom cell 
representation, described above, for interfaces which have contiguous mesh points. 
The disadvantages of this approach are that the interpolation is essentially linear, 
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Non-Contiguous Mesh Block Interface 
Boundary Coupling Scheme 


Mesh Block Interface Structure 



© 


Mesh Block # 1 Grid Lines 
Mesh Block #2 Grid Lines 
Mesh Block #i Phantom Cell of Interest 

Mesh Black *1 Phantom Cell . Face Cent" 

Mesh Block #2 Boundary Cell Face Cente 

Electrical Circuit Analogy 
A 

Resistance Element Proportional to Spacing 


B 



Figure 3 


.10: ADPAC-AOACR Non-Contiguous Mesh Block Coupling Scheme 
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and that the scheme is not uniformly conservative, and is thus subject to errors is 
mass, momentum, and energy conservation. These sources of error scale roughly with 

the mesh spacing, and therefore, the finer the mesh , the lower the impact these errors 
have on the overall solution. 

The third type of mesh block interface boundary patching scheme relates to 
steady flow calculations involving turbomachinery components with multiple blade 
rows, or to a lesser extent, steady flow calculations involving mesh systems with differ- 
ing circumferential extents. This approach is a technique for the prediction of steady 
flow through multistage turbomachinery based on the mixing plane concept. A mixing 
plane is an arbitrarily imposed boundary inserted between adjacent blade rows across 
which the flow is “mixed out” circumferentially. This circumferential mixing approx- 
imates the time-averaged condition at the mixing plane and allows the aerodynamic 
solution for each blade passage to be performed in a steady flow environment. The 
mixing plane concept was recently applied to realistic turbofan engine configurations 
by Dawes [6]. Flow variables on either side of the mixing plane are circumferentially 
averaged and passed to the neighboring blade row as a means of smearing out the 
circumferential nonuniformities resulting from dissimilar blade counts. The mixing 
plane concept is a cost-effective approach computationally because the flow is steady, 
and the individual blade passage domains are limited to a near-blade region. Unfortu- 
nately, the accuracy of this approach is clearly questionable under some circumstances 
because of the placement of the mixing plane and the loss of spatial information re- 
sulting from the circumferential averaging operator. Numerically, then, the mesh 
points in two adjacent blocks at a mixing plane interface must have identical axial 
and radial coordinates. This permits specifying the circumferential distribution of 
phantom cells in one block as the circumferential average of the near boundary cell 
data from the adjacent block. This concept is illustrated graphically in Figure 3.11. 
The individual primitive flow variables are each circumferentially averaged, and used 
directly in the phantom cell data specification. Attempts to use several of the more 
complex mixing analyses based on the circumferentially averaged flow data, or the 
use of characteristic-based boundary specification based on the flow data were both 
subject to numerical instabilities under particular circumstances, and therefore, the 
direct specification method is used instead. 


45 


Mixin g Plane Concept Mesh Block 
Boundary Coupling Scheme 


Geometry and Mesh Block Structure 



Mesh Block # 1 Grid Lines 
Mesh Block #2 Grid Lines 


Mixing Plane Interface 


Mesh Block #1 


Mesh Block #2 




Circumferential average of near-boundary 
data in Mesh Block #2 is applied uniformly 

across the circumferential distribution of 
it > / „ ^ if 7 


Figure 3.11: 


ADPAC-AOACR Circumferential Averaging (Mixing Plane) Mesh 
Block Coupling Scheme for Multiple Blade Row Calculations 
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he fourth type of mesh block interface boundary patching scheme relates to 
time-dependent calculations involving relative motion between adjacent mesh blocks. 
In this case, the mesh interface is essentially noncontiguous, and the relative posi- 
tion between mesh cells is changing in time. This type of boundary scheme is used 
or the prediction of rotor/stator interaction for multiple blade row turbomachines, 
his type of interface is complicated by the fact that the interface may be defined 
y different mesh blocks at different times based on the degree of motion between 
the adjacent mesh blocks. For example, if a calculation for a single stage compressor 
involves 3 passages of the upstream blade row, and four passages of the downstream 
blade row, then it is obvious that at some point during the calculation, each of the 
t ree mesh blocks in the upstream blade row representation will at some time have 
to communicate with each of the four mesh blocks in the downstream blade row rep- 
resentation. The specific mesh block in the communication path is determined by the 
degree of rotation in both the upstream and downstream blade rows, and is inher- 
ently time dependent. This task is simplified considerably by requiring that the mesh 
points in each block along the interface have constant axial and radial coordinates. 
This simplification reduces the overall problem to a time-space interpolation in the 
circumferential direction only (see Figure 3.12). 

With this simplification, it is possible to simply define the interpolation data by 
resolvmg the relative position between the mesh blocks, extracting the near bound- 
ary data values from the neighboring mesh blocks, and interpolating for the phantom 
cell data value based on the current circumferential position of the phantom cell face 
center relative to the neighboring interpolant data arrays. Without modification, this 
scheme is nonconservative and can result in spurious errors in mass, momentum, and 
energy conservation across the mesh block interface. In order to rectify this problem, 
a modified interpolation scheme was tested which attempted to enforce global conser- 
vation of the interpolated data across each circumferential interpolation “slice” of the 
mesh. This restriction was enforced by globally summing the interpolated values for 
each circumferential interpolation and modifying the summation to satisfy the global 
conservation defined by the interpolant data. The individual interpolated phantom 
cell values were then uniformly modified to reflect the global conservation correction. 
Numerical experiments with this scheme compared to the unmodified scheme were 
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Time/Space Resolved Interpolation Mesh 
Block Boundary Coupling Scheme 

Geometry and Mesh Block Structure 

(Circumferential Pitch of Each Blade Row 
Representation is Constant) 



Blade Row #1 


Boundary Surface Interface 

Blade Row #2 




distribution of phantom cells in Blad 

Figure 3.12: ADPAC-AOACR Time/Space Solved Mesh Block Coupli„ g Scheme 

for Multiple Blade Row Calculations 




virtually identical, and it was not possible to draw any hard conclusions regarding the 
advantages of one scheme over another. Following this test, the correction scheme was 
eliminated for simplicity. The effects of conservation across block interface boundaries 
is likely to be a debatable issue; however, there is considerable numerical evidence 
available which suggest that solution accuracy (for the types of calculations consid- 
ered in this study) is not significantly influenced by non-conservative interpolation 
schemes ( [12], [39]) and that supposed conservative interface coupling schemes [13] 
do not necessarily ensure conservation when used with numerical schemes which are 
not conservative in time. 


3.13 Solution Procedure 

The overall solution procedure begins by defining a set of initial data, and ad- 
vancing the solution from that point forward in time until the desired solution (steady 
state, time-periodic, or finite time interval) has been reached. Initial data is normally 
specified as a uniform flow, or may be read in as a “restart” of a previous existing so- 
lution. Normally, for steady flow calculations, the “full” multigrid startup procedure 
is utilized to accelerate convergence by initializing the solution on a coarse mesh be- 
fore incurring the expense of fine mesh iterations. Steady state solutions are deemed 
converged when the average residual R has been reduced by a factor of 10 - ^, or 
when the residual has ceased to be reduced. It is possible that for some steady flow 
calculations, the solution is truly unsteady (i.e. - vortex shedding behind a circular 
cylinder) and in these cases the residual may not be reduced beyond a certain limit. 

For time-dependent flow simulations, it has been found that it is often desirable 
to first initialize the solution in a steady flow manner using mixing planes or other 
time-averaged boundary condition procedures to define a reasonable approximation to 
the time-averaged flow before proceeding with the time-accurate iteration procedure. 
The unsteady solution may then be advanced in time with a uniformly specified time 
step, until a time-periodic solution is achieved. The amount of time required to reach 
a time-periodic solution is typically highly case dependent, and such solutions must 
be monitored closely to determine when the solution is satisfactory. 

Some comments are in order to explain the solution procedure for high speed 
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rotors, many of which are iUustrated in the results section in this report. It has been 
found that solutions for high speed compressor rotors typically require exit static 
pressure in excess of the inlet total pressure. As a result, it is often not possible to 
obtain an immediate solution under the desired operating condition because of the 
potential for reversed flow at an exit boundary which results from the relatively large 
exit static pressure. To overcome this problem, the following procedure has been 
implemented to obtain high speed rotor flowfield solutions. Typically, the solution is 
started using an exit static pressure which is lower than the inlet total pressure. For 
example, for a fan rotor with a design exit static to inlet total pressure ratio of 1.13, 
the solution may be started with an exit static to inlet total pressure ratio of 0.95. The 
resulting solution typically chokes the flow in the blade passage and sets up a shock 
system which is consistent with the choked flowfield. Once this flow is established, 
the exit static pressure may be raised incrementally, and the solution restarted several 
times until the desired exit static pressure ratio is reached. It should be emphasized 
that the increases in back pressure must be done conservatively to avoid prematurely 

stalling the blade passage. 
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4 . RESULTS 


In this chapter, numerical results from the ADPAC-AOACR 3D Euler/Navier- 
Stokes code are presented and compared with existing data to demonstrate and verify 
the accuracy of the analysis. The results are primarily based on applications involv- 
ing ducted and unducted fans with multiple blade rows; however, several additional 
calculations are presented to demonstrate the versatility of the code. 

Steady flow predictions involving turbomachinery geometries are presented first, 
followed by several calculations of time-dependent flowfields involving multiple blade 
row turbomachines. 


4,1 2-D Compressor Cavity Flow (Steady Flow) 

In light of the recent addition of the two-dimensional block solution capability 
in the ADPAC-AOACR analysis code, a two-dimensional test demonstration for a 
complex flow is appropriate. One such flow of interest to turbomachinery designers 
is the flow in a compressor hub cavity for compressors which employ inner banded 
stators. The exact characteristics of the leakage flow through the cavity and the 
resulting interplay between the leakage flow and the primary gas flow are not com- 
pletely understood. As such, inner banded stator seal design is currently an area of 
concentrated investigation for high performance gas turbine engines. 

The ADPAC-AOACR analysis was applied to predict the two-dimensional flow 
in an inner-banded stator hub cavity design considered for an advanced compressor 
design at Allison. This geometry, and the 149x61 mesh system used in the analysis 
are illustrated in Figure 4.1. This complex mesh system is actually rectangular in 
computational space, and the application of the appropriate boundary conditions is a 
trivial matter for the ADPAC-AOACR analysis. The time-marching calculation was 


51 




Figure 4.1: Inner Banded Stator Cavity Geometry and ADPAC-AOACR 2-D Mesh 

System 


advanced with a CFL number of 5, utilized 3 levels of multigrid, and employed the 
full multigrid startup procedure. The residual convergence history resulting from this 
calculation is given in Figure 4.2. The full multigrid startup procedure consisted of 
200 iterations on each of the 2 coarse mesh levels, foUowed by 100 fine mesh iterat.ons. 
The overall convergence for this case was excellent, with the average residual reduced 
a full four orders of magnitude in only 100 fine mesh iterations. 

The predicted steady flow velocity vector pattern for this flow is given in Fig- 
ure 4.3. It should be noted that the flow in the cavity is from left to right in this 
orientation. The flow consists of a complex pattern of recirculating vortices, coupled 
with the high speed, high loss flow through the knife seals. It is particularly interesting 
to note the standing vortices at the inlet and exit of the cavity flowpath. Although 
no experimental data were available to verify these results, the predicted 2-D flow 
was in good agreement with a similar 3-D flow prediction from the ADPAC-AOACR 
analysis, and a similar 2-D flow predictions based on a separate 2-D Navier-Stokes 

code. 
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Banded Stator Cavity Flow Calculation 
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Figure 4.3: ADPAC-AOACR Predicted Velocity Vectors for 2-D Inner Banded Sta- 

tor Cavity Flow Analysis 


4.2 Model Counterrotating Propfan (Steady Flow) 

The first results presented are for a model counterrotating propfan geometry. 
Geometric and aerodynamic design parameters for the model counterrotat.ng pro- 
peller are given in Figure 4.4. The calculations presented here were based on a cru.se 
condition with a flight Mach number of 0.72, and advance ratio of 3.14. The forward 
and rearward blade setting angles were 57.4 and 60.0 degrees, respect.vely. Steady 
flow performance predictions were obtained for this case by solving the Naver-Stokes 
equations with the inter-blade row mixing plane concept described m Section . . 
The grid therefore consists of two blocks (one for each blade row), divi eye 
mixing plane approximately midway between the two blade rows. Our experience sug- 
gests that propfan airfoils are sufficiently flexible that in order to accurately capture 
aerodynamic performance, it is necessary to perform coupled structural/aerodynam.c 
analyses to deduce the airfoil shape under the desired operating conditions For e 
purposes of this demonstration, the individual mesh blocks were generated based on 

an assumed fixed running airfoil shape. 
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It should be noted that this grid system corresponds to standard configuration 
#7 as described in the ADPAC-AOACR Computer Program Users Manual [9j. The 
grid block sizes were 65x49x33 and 57x49x33, respectively. An illustration of the 
mesh system is given in Figure 4.5. The calculations presented here are essentially 
for demonstration purposes only. In practice, the exit boundary should probably be 
extended farther downstream to obtain best results. 

Calculations for this geometry utilized 3 levels of multigrid, and employed the 
full multigrid startup procedure. A total of 100 iterations were performed on each of 
the coarse mesh levels during the full multigrid startup, followed by 200 iterations of 
the fine mesh solver. The CFL number was fixed at 5.0, and the modified four stage 
time-marching algorithm was employed. An illustration of the convergence history 
for this case is given in Figure 4.6. 

The overall convergence of the multigrid procedure was observed to be quite 
rapid, requiring only 200 fine mesh iterations to achieve a 3 order reduction of the av- 
erage solution residual (error). An illustration of the predicted surface static pressure 
contours and radial blade-to-blade static pressure contours are given in Figure 4.7. 
The blade to blade static pressure contours clearly illustrate the tip vortex flow near 
the suction side of each blade. No leading edge vortex was observed in the predictions 
for either blade row under this operating condition, which appears to be supported 
by experimental evidence. The overall predicted power coefficient for this case was 
2.05 which is in good agreement with the design value of 2.035. 

Two images are presented in Figures 4.8 and 4.9 to illustrate the effectiveness 
of the mixing plane boundary patching procedure for this type of flow. Figure 4.8 
illustrates predicted absolute Mach number contours for the axisymmetnc averaged 
flowfield. The blade boundaries and location of the mixing plane are outlined. The 
corresponding axisymmetric averaged flow static pressure contours are given in Fig- 
ure 4.9. Both sets of contours display a smooth transition across the mixing plane 
and display no noticeable abnormalities as a result of the mixing plane assumption. 
Another plot, given on Figure 4.10, illustrates the circumferential blade to blade 
static pressure contours at midspan for the mixing plane calculation. The mesh block 
boundaries are outlined, and the mixing plane boundary is visible between the blade 
rows. Naturally, the contours do not mate at the mixing plane boundary because 
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Mesh Block Structure: 

1 Mesh Block/Blade Row (2 Mesh Blocks Total) 
Mixing Plane Approximately Midway Between 
Blade Rows 

Number of Blades in Forward Row: 5 

Number of Blades in Rearward Row: 5 


Design Flight Mach Number: 0.72 

Design Power Coefficient 2.035 

Design Advance Ratio: 2.94 

Design Tip Speed: 750 ft/s 

Power Loading/D**2: 37.15 



Figure 4.4: Model Counterrotating Propfan Model Geometric and Aerodynamic De- 

sign Parameters 
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Figure 4.5: ADPAC-AOACR Two-Block 

Model Counterrotating Propfai 
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Log 10 (Residual) 


Model Counterrotating Propfan 



Figure 4.6: ADPAC-AOACR Multigrid Convergence History for Model Counterro- 

tating Propfan Steady Flow Analysis Based on Mixing Plane Concept 
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Figure 4.7: ADPAC-AOACR Predicted Steady Surface Static Pressure Countours 

and Post Rotor Static Pressure Contours for Model Counterrotating 

Propfan 
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of the circumferential averaging procedure. It is apparent, however, that the mixing 
plane does not unduly distort the flowfield contours in the vicinty of either blade row. 

This calculation clearly demonstrates one application which can benefit from 
the accuracy and cost-effectiveness of the mixing plane approach. Future studies 
concentrating on defining the limits (in terms of interblade row spacing) would be 
helpful to define the geometric regime where the mixing plane assumption is valid. 


4.3 NASA 1.15 Pressure Ratio Fan (Steady Flow) 

In this section, calculations are presented for the viscous flow through a 1.15 
pressure ratio fan stage originally tested by NASA [42]- [47] and utilized extensively 
under this contract for analysis [21], [7]. Descriptions of the geometry and design 
parameters for the NASA 1.15 pressure ratio fan are given in Fig. 4.11. 

The NASA 1.15 pressure ratio fan is representative of a 25:1 bypass ratio tur- 
bofan engine fan stage, and therefore closely approximates the ducted propfan con- 
cept propulsion system. Several steady flow calculations were performed to permit a 
comparison of predicted results with the high speed experimental data published in 
Ref. [42]. 

A meridional view of the geometry and steady flow grid system are given in 
Fig. 4.12. 

The steady flow calculations for this geometry were performed using the mixing 
plane inter-blade row model between the fan rotor and stator. It should be noted that 
this grid system corresponds to standard configuration $8 described in the ADPAC- 
AOACR Computer Program Users Manual (see [9]). The mesh system, illustrated in 
Figure 4.12, consists of 4 mesh blocks, two each associated with the rotor and the 
stator, respectively. The overall mesh block sizes were 89x37x33, 89x25x33, 81x37x33, 
and 81x25x33. The rotor was represented by 33 points in the radial direction, result- 
ing in 5 points to define the clearance region above the rotor tip. Calculations are 
presented here for a flight Mach number of 0.75, and a fan rotational speed of 9167 
rpm (100% speed). The calculation utilized 3 levels of multigrid and employed the 
full multigrid startup procedure. A total of 100 iterations were performed on each of 
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Figure 4.8: ADPAC-AOACR Predicted Steady Axisymmetric Averaged Mach Num- 

ber Contours for Model Counterrotating Propfan 
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Figure 4.9: ADPAC-AOACR Predicted Steady Axisymmetric Averaged Static Pres 

sure Contours for Model Counterrotating Propfan 








(dimensions in cm) 


Figure 4.11: NASA 1.15 Pressure Ratio Fan Geometry and Design Parameters 
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Figure 4.12: ADPAC-AOACR Four Block Mesh System for NASA 1.15 Pressure 

Ratio Fan 
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the coarse mesh levels during the full multigrid startup, followed by 300 iterations of 
the fine mesh solver. 

The overall convergence of the multigrid procedure was very good for this calcu- 
lation, in spite of the added complexity of the multiple block discretization and the 
mixing plane coupling scheme. Predicted axisymmetric averaged flow contours for 
absolute Mach number and static pressure are illustrated on Figures 4.13 and 4.14, 
respectively. 

The cowl, spinner, blades, and mixing planes are outlined on both figures. The 
axisymmetric flow displays a smooth transition across the mixing plane, although 
this feature was certainly not unexpected considering the axial separation between 
the blade rows for this case. 

A comparison of the predicted stator exit radial total pressure distribution with 
experimental data and numerical data obtained with the average-passage formulation 
developed under Task I of this contract [21] is presented in Figure 4.15. The inviscid 
prediction using the HPR03D developed under Task 1 clearly overpredicts the total 
pressure ratio at the nozzle entrance as a result of the lack of viscous losses. The 
present prediction shows good agreement across the radial span excepet near the 
hub, where the ADPAC-AOACR results clearly underpredict the total pressure rise 
through the machine. This discrepency is thought to be due to numerically generated 
losses as a result of poor grid refinement near the spinner leading edge. These losses 
are ultimately convected downstream and end up showing up along the hub at the 
nozzle entrance. An illustration of the three-dimensional flowfield for this case is 
given in Figure 4.16. This figure displays predicted surface static pressure contours 
for a 3-D representation of the entire fan geometry, with a portion of the fan cowl cut 
away to expose the inner fan detail. 

Additional calculations for this geometry were performed under various flow con- 
ditions (although the results will not be reported in detail here) including a supersonic 
flight test case at a freestream Mach number of 2.5. (In this case, the supersonic flow 
condition was generated by specifying an initial freestream Mach number of 2.5, and 
simply neglecting the application of any boundary condition at the mesh inlet, im- 
plying that the initial values are never modified, which provides a physically realistic 
supersonic inflow boundary treatment.) This case is mentioned only as an example 
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Figure 4.13: ADPAC-AOACR Predicted Axisymmetric Averaged Flow Absolute 

Mach Number Contours for NASA 1.15 Presure Ratio Fan 
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Figure 4.14: ADPAC-AOACR Predicted Axisymmetric Averaged Flow Static Pres 

sure Contours for NASA 1.15 Presure Ratio Fan 
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Figure 4.16: ADPAC-AOACR Predicted 3-D Surface Static Pressure Contours for 

NASA 1.15 Pressure Ratio Fan (M=0.75) 
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of the flexibility of the code for solving a wide variety of flow problems. In most cases 
tested, the multigrid procedure rapidly converged to a steady state solution (3 order 
of magnitude reduction in the average fine mesh residual) in less than 300 fine mesh 

iterations. 


4.4 GMA3007 Fan Section (Steady Flow) 

Steady flow verification of the ADPAC-AOACR code was also performed for an 
advanced turbofan engine fan rotor geometry taken from the Allison GMA3007 tur- 
bofan engine. The Allison GMA3007 is a 5:1 bypass ratio turbofan engine which 
produces 7,000 pounds of thrust. The fan section geometry considered consists of 
three blade rows including the fan rotor, bypass stator, and core stator. The flow 
is complicated by the presence of the core flow splitter, and the fact that there are 
multiple exits (bypass flow, core flow) for which proper boundary data must be spec- 
ified. An outline of the GMA3007 fan section geometry is given in Figure 4.17. This 
domain was arbitrarily decomposed into three mesh blocks, as shown in Figure 4.17. 

The individual mesh blocks correspond to one of the three blade rows in the fan 
section. For steady flow analysis, the mixing plane concept described in Section 3.2 
was invoked, and therefore the circumferential pitch of each mesh block corresponds 
to the periodic spacing of the included blade row. In addition to the axial placement 
of the mixing planes between blade rows, the mesh system was designed to illustrate 
a radial mixing plane which separates the mesh blocks for the bypass stator (#2) and 
the core stator (#3) just upstream of the core flow splitter. Obviously, proper place- 
ment of the mesh block boundaries could have removed the need for this interface, 
but a demonstration of the radial mixing plane was thought to be justified in light 
of future tasks under the current contract which wiU utilize this type of boundary m 

more detail. 

The actual mesh system is detailed in Figures 4.18 and 4.19. The grid block sizes 
were 85x73x29, 101x33x29, and 101x41x29, respectively. The grid was constructed 
by using the program TIGG3D , and special care was taken to maintain a reasonable 
clustering of mesh points near solid boundaries to resolve the viscous shear layers, and 
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Mesh Block Structure 



Figure 4.17: Allison GMA3007 Fan Section Geometry 
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Figure 4.18: Axisymmetric Detail of ADPAC-AOACR Three-Block Mesh System for 

Allison GMA3007 Fan Section 


to maintain grid indices which were compatible with the multigrid solution technique. 
The fan rotor includes a tip clearance region which was discretized with 9 points 
radially to resolve the tip clearance flow. 

The numerical solution was performed using 3 levels for multigrid and was ini- 
tiated with the full multigrid startup procedure. As is typical for high speed rotor 
flowfield calculations, the solution must be initiated at a back pressure which is lower 
than the intended solution back pressure. The usual procedure is to slowly raise the 
back pressure to the desired setting by restarting the calculation a number of times 
from a previous run. This procedure prevents excessive shock motion while setting 
up the flowfield which can prematurely cause the flowfield to stall. In this case, since 
there are multiple domain exit regions, following the recommended startup procedure 
is even more important to eliminate large oscillations in the flow split between the 

core and bypass ducts. 
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Figure 4.19: Circumferential Detail of ADPAC-AOACR Three-Block Mesh System 

for Allison GMA3007 Fan Section 


Following the startup procedure described above, a design point (100% speed) 
solution was achieved in approximately 1500 fine mesh iterations. Experimented data 
were available at several points throughout this flow. 

Figure 4.20 illustrates a comparison of the predicted and experimental fan ro- 
tor exit radial total pressure ratio distribution taken approximately one half chord 
downstream of the rotor. Predicted results from several other calculations for this 
geometry are also included on this Figure to illustrate the accuracy of the ADPAC- 
AOACR code compared with other analyses. It should be mentioned that some slight 
differences in operating condition (mass flow rate, etc.) existed between the various 
predictions displayed on this figure, which is at least partially responsible for the 
differences in the results. The ADPAC-AOACR results are clearly in good agreement 
with the experimental data at this point, and is at least as accurate as any other 
analysis which has been tested for this geometry to date. An interesting observation 
from this result is that the ADPAC-AOACR prediction is in good agreement with 
the ADPAC-APES code, which was also developed under this contract, and which 
utilizes the more complicated average-passage solution scheme for multiple blade row 
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Figure 4.23: Predicted Surface Static Pressure Contours for Allison GMA3007 Fan 

Section 


78 


COLOR PHOTOGRAPH 



turbomachines. One would expect that the average-passage equation system provides 
a more accurate model for steady flow analysis of interblade row interactions than 
the simpler mixing plane concept; however, the large axial spacing between blade 
rows in this turbomachinery geometry suggests that this advantage is likely to be ess 
influential, and the good agreement between predictions supports this observation. A 
similar comparison of the predicted and experimental bypass stator leading edge ra- 
dial total pressure and total temperature distributions is given on Figure 4.21. Once 
again, the ADPAC-AOACR results are in good agreement with the average passage 

prediction for this geometry. 

An illustration of the predicted axisymmetric-averaged static pressure contours 
for this calculation are presented in Figure 4.22. The stagnation region around the 
core flow splitter is clearly defined, and it is clear that there is some slight negative 
incidence in the vicinity of the splitter leading edge. This incidence is clearly depen- 
dent on the exit pressure chosen for both the core flow and bypass flow domains, and 
our experience has been that the solution is relatively sensitive to small changes in 

either specification. 

Finally, to illustrate the three-dimensional nature of this flowfield, an illustra- 
tion of the predicted static pressure contours for a complete representation of the 
GMA3007 fan section is given in Figure 4.23. The contours clearly display the rotor 
shock location, and the flowfields associated with the core and bypass stators. The 
complexity of this figure demonstrates the capability of the ADPAC-AOACR code to 
analyze complex, realistic turbofan engine flowfields. 

4.5 Oscillating Flat Plate Cascade (Unsteady Flow) 

Before proceeding with a discussion of time-dependent predictions for multi- 
ple blade row turbomachines, a demonstration of the time-accuracy of the ADPAC- 
AOACR solution procedure seems appropriate. A numerical stud} was performed 
involving the time-dependent flowfield about an oscillating flat plate cascade. This 
calculation serves to further verify the time accuracy of the ADPAC-AOACR code, 
and can provide some insight into the possibility of utilizing nonlinear time-marching 
numerical solution techniques for the prediction of blade motion induced acoustic 
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phenomena in turbomachinery components. A detailed discussion of calculations of 
this type was recently presented by Huff [45], and these results were used to guide 
the present calculations. 

The calculation was based on a zero stagger flat plate cascade with a pitch to 
chord ratio of 1.0. The cascade motion was represented by a pitchwise plunging 
harmonic oscillation. This configuration is illustrated in Figure 4.24. The cascade 
motion was programmed into the ADPAC-AOACR code by temporarily incorporating 
a time- varying blade rotation term. The solution was initiated from a steady, uniform 
flowfield, and advanced in time under the influence of the harmonic blade motion, 
until a time-periodic solution was achieved. Calculations were performed on both a 
fine and coarse mesh system for a reduced frequency of 4.0 (defined as 


k = 


UJC 

Wi) 


(4.1) 


where k is the reduced frequency, w is the harmonic oscillation frequency, c is the 
flat plate chord, and V i is the freestream velocity. The freestream flow Mach number 
was 0.5. The amplitude of the cascade vibration was 0.1% of the blade chord. The 
mesh systems were actually based on a 2-D grid system generated by personnel at 
the NASA-Lewis Research Center, which was stacked radially to construct a 3-D 
mesh system. The 3-D mesh system was constructed with a hub to tip ratio of 
0.996 and a radius to chord ratio of 100.0, which were believed to be sufficient to 
minimize radial flow gradients and non-planar flow features resulting from the use of 
a cylindrical coordinate system.. The coarse grid was decomposed from the fine grid 
by eliminating every other grid line. An illustration of the blade-to-blade component 
of both mesh systems is given in Figure 4.25. 

The fine mesh size was 121x5x41, while the coarse mesh was 61x5x21 (the radial 
distribution of mesh nodes is kept low (5) since the flow is essentially two-dimensional 
(no radial flow gradients)). The calculations were advanced over approximately 25 
cycles of the plunging cascade oscillation using a constant CFL number of 1.0. The 
calculation was performed on an IBM RS6000 Model 540 workstation and required 
approximately 28 hours (estimated) to complete the fine mesh solution. The time- 
periodicity of the solution at this point was verified by comparing local static pressure 
time-histories at several points along the cascade surface and verifying that the local 
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Figure 4.24: Oscillating Flat Plate Cascade Geometry 
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GEOMETRY 


Flat Plate Cascade Harmonic Oscillation Study 


Coarse Mesh ( 61x5x21 --> axial, radial, circumferential ) 



Fine Mesh (121x5x41 --> axial, radial, circumferential ) 
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Figure 4.25: Oscillating Flat Plate Cascade Blade to Blade Mesh Systems 
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pressure was periodic with respect to the known oscillation frequency. During the 
calculation, time histories of the local circumferential static pressure differential across 
the flat plate airfoils were stored. This static pressure differential, defined as 

Ap(z,r) = p(x,r,6) — p(x,r,0 + AO) (4-2) 


where A 0 is the cascade pitch, is representative of the instantaneous local aerody- 
namic loading on the plate surface. The final cycle of these time histories was decom- 
posed into a Fourier series, and the first harmonic component of the Fourier series 
was then plotted as real and imaginary components of the plate pressure differential 
response to the cascade motion. The predicted real and imaginary components of the 
flat plate pressure differential response are illustrated in Figures 4.26and 4.27. Note 
that the response is illustrated in terms of a pressure coefficient as: 


A C p = 


Ap 

pv^h k 


(4.3) 


where h is the amplitude of oscillation. 

The predicted results were compared with results from the linear theory presented 
by Smith [44]. The comparison of the real and imaginary components for the plate 
pressure differential response given in Figures 4.26 and 4.27 includes both coarse and 
fine mesh ADPAC-AOACR Euler predictions and the Smith results. This comparison 
clearly suggests reasonable agreement with the linear theory, except in the vicinity 
of the leading edge, where nonlinear effects could be more important, in which case 
the linear theory is likely to be less accurate. Due to a lack of time, no specific 
measures were taken (such as reducing the amplitude of oscillation) to evaluate the 
magnitude of this apparent nonlinear behavior. The predicted results demonstrated 
a marked improvement in the correlation with the Smith results with increased grid 
resolution. One benefit of the present mesh system is that the results are likely to 
be fairly insensitive to inlet or exit boundary conditions It should be observed that 
the present calculations are likely to be relatively insensitive to the inflow/outflow 
boundary conditions as a result of the coarsening of the axial mesh spacing away 
from the flat plate cascade. This coarsening effectively damps out the unsteady waves 
in the farfield and minimizes nonphysical reflections from the far field boundaries. 
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Figure 4.26: Comparison of y4DP.4C-.40/lCi? Prediction and Smith Linear Theory 

for Real Component of Airfoil Surface Pressure Response for the Os- 
cillating Flat Plate Cascade Test Case 
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Figure 4.27: Comparison of ADPAC-AOACR Prediction and Smith Linear Theory 

for Imaginary Component of Airfoil Surface Pressure Response for the 
Oscillating Flat Plate Cascgde Test Case 
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Additional details on the effects of mesh coarsening and inflow/outflow boundary 
treatments are detailed in Ref. [45]. 

It is recommended that additional calculations of this type be performed in the 
future for more complicated flow cases, and mesh systems involving multiple blocks 
to analyze the influence of the interblock communication boundary conditions on the 
unsteady flow results. 


4.6 SR7 Unducted Propfan - Cylindrical Post Interaction Study 

(Unsteady Flow) 

In this section, an application of the ADPAC-AOACR program to predict the 
time-dependent flowfield resulting from the aerodynamic interaction between a sta- 
tionary cylindrical post and a rotating unducted propfan is described. These results 
represent the first presentation of the time-dependent multiple blade row flow anal- 
ysis capability of the ADPAC-AOACR program. The flow conditions for this case 
were based on results from the experimental study described by Bushnell et al. [40]. 
Bushnell presents wind tunnel measurements of a single rotation propfan mounted 
behind a stationary cylindrical post. The purpose behind this study was to deter- 
mine the effect of the interaction between the wake from the cylindrical post and the 
rotating propfan airfoil. The experimental configuration is shown diagramatically in 
Figure 4.28. The propfan design utilized in this study was based on the SR7 propfan. 
Geometric and aerodynamic design parameters for the original SR7 propfan are given 
in Figure 4.29. 

The experimented measurements included steady state airfoil static pressure data, 
and time- dependent airfoil static pressure data resulting from both angle of attack, 
and the unsteady interaction with the cylindrical post. The steady state and experi- 
mental pressure data locations are given in Figure 4.30. 

Because of wind tunnel power limitations, the number of propfan blades used in 
the experimented study was reduced from 8 (the original SR7 design blade count) to 2. 
This reduction in blade count resulted in a fortuitous simplification in the numerical 
analysis of the unsteady interaction between the post and the rotating propfan. Since 
the effective blade count for the propfan and the cylindrical post are identical (in this 
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SR7 Propfan/Cylindrical Post Aerodynamic 
Interaction Study Configuration 


Stationary 



Figure 4.28: SR" Preplan - Cylindrical Port Interaction Study Experimental Con- 

figuration 
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Characteristic 

SR-7A 

Number of blades 

8 

lip sweep angle,* deg 

41 

Model diameter, cm (In.) 

62.23 (24.5) 

Tip speed, m/sec (ft/sec) I at 10.669 m(35 000 ft) 

243.8 (800) 

Power loading, kW/m? (shp/ft 2 ) jl.S.A. altitude 

256.05 (32) 

Activity factor, AF 

227 

Integrated design lift coefficient,^ 

0.^02 

A1 rf ol 1 s 

NASA 16 and 65/CA 

Ratio of nacelle maximum diameter to propeller diameter 

0.35 

Cruise design Mach number 

0.80 

Cruise design advance ratio 

3.06 

Cruise design power coefficient 

!_!!_ 


a Geometr1c measurement from planform. 


Figure 4.29: Geometric and Aerodynamic Design Parameters for the Original SR7 

Propfan 
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case the 

stationary cylindrical post is modeled as a blade row with 2 blades), the resulting 
mesh system requires only one mesh block per blade row. An illustration of the mesh 
block structure for this case is given in Figure 4.31. It should be mentioned that tbs 
mesh structure corresponds to Standard Configuration #7 described in the ADPAC- 
AOACR Computer Program Users Manual [9]. 

During the course of this numerical study, several mesh systems were tested 
with the ADPAC- AOACR code to predict the time-dependent static pressure on the 
propfan airfoil as it interacts with the wake from the cylindrical post. Grid resolution 
was immediately identified as a problem with most mesh systems tested. In order 
to adequately convect the wake from the cylinder, a large number of circumferential 
mesh points were required in both grid systems. This was predominantly evident m 
the propfan mesh block. As the wake is convected into the downstream mesh block, 
the circumferential spacing of the mesh typically varies dramatically based on the 
relative position of the profan blades as they rotate during the time-accurate solution 
process. Due to computational costs, it was not feasible to maintain the required 
circumferential spacing in the propfan mesh block required to accurately convect the 
wake throughout the blade rotation. The best compromise appeared to be to move 
the interface between the two mesh blocks as close to the propfan as possible, and to 
allow the stationary mesh block to convect the wake as far downstream as possible. 
This suggests that the influence of the wake on the propfan airfoil is only important 
when the airfoil and the wake are in close proximity, as the wake is effectively lost in 
the propfan mesh when the airfoil is rotated away from the plane of the wake. 

An axisymmetric view of the final mesh system used to provide the numerical 
results is given in Figure 4.31. This mesh system provided a reasonable compro- 
mise between solution accuracy, and computational cost for the prediction of this 
complicated unsteady flow. 

The time-dependent calculation was initiated from a steady state solution on 
the same grid using the mixing plane concept described in Section 3.2 to couple the 
adjacent, relatively rotating blade rows. Once the steady flow solution was reached, 
the time-dependent calculation was initiated using the space/time resolved interblade 
row coupling procedure described in Section 3.3. The solution was advanced in time 
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Figure 4.30: SR7 Propfan - Cylindrical Post Aerodynamic Interaction Experimental 

Airfoil Surface Data Locations 






for 3 complete revolutions of the propfan rotor. It was observed that for all practical 
purposes, a time-periodic solution had been achieved after the first revolution. This 
is most likely a result of the fact that the cylinder and propfan airfoil count is so low, 
and the airfoils actually behave as if they are isolated, rather than cascaded. The 
time-periodic results presented here pertain to the final revolution of the rotor. 

Predicted propfan surface pressure histories given as the unsteady deviation from 
the steady state static pressure divided by the steady state static pressure are com- 
pared with time-resolved experimental data in Figures 4.32-4.35. 

The predicted aerodynamic response is in good agreement with the experimen- 
tal traces over all of the measurement locations surveyed. In most cases, the peak 
pressure differential resulting from the interaction between the propfan and the wake 
is slightly underpredicted by the ADPAC-AOACR code. It is also evident from the 
predicted time-interval over which the airfoil static pressure history is disturbed that 
the wake width has been overpredicted. It is believed that both of these effects could 
be improved by increasing the overall mesh resolution (although at the expense of 

increasing the overall computational cost). 

An illustration of the predicted instantaneous total pressure contours at the 
midspan of the propfan airfoil is given in Figure 4.36. This figure clearly depicts 
the propfan airfoil slicing through the cylinder wake, and the resulting effects on the 
propfan airfoil pressure field. An interesting side note for this case is that the ex- 
perimental Reynolds number for the cylinder was within the range (50 < Re D < 
400000) where Karman vortex streets occur. Graphical visualization of the instanta- 
neous flow behind the cylinder at several times did not depict any vortex shedding 
from the cylinder. It is unlikely that the numerical grid used in this calculation is 
truly dense enough to accurately support any physical vortex shedding phenomena. 


4.7 Model Counterrotating Propfan (Unsteady Flow) 

In this section, numerical results are presented for a time-dependent calculation 
of the flow through the model counter rotation propeller described in Section 4.2. 
Steady flow predictions for this geometry based on the mixing plane concept were 
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Figure 4.32: 


Comparison of Time-Resolved Predicted and Experimental Propfan 
Surface Static Pressure Histories for SR7 Propfan - Cylindrical Post 
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Figure 4.33: Comparison of Time-Resolved Predicted and Experimental Propfan 

Surface Static Pressure Histories for SR7 Propfan - Cylindrical Post 
Aerodynamic Interaction Study (36.7% Axial Chord, 64.1% Radial 
Span, Suction Side) 94 
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Figure 4.34: Comparison of Time-Resolved Predicted and Experimental Propfan 

6 Surface Static Pressure Histories for SR7 Propfan - CylindncaJ Pos 

Aerodynamic Interaction Study (10.0% Axial Chord, 64.1% Radial 
Span, Pressure Side) 95 
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Figure 4.35: Comparison of Time- Resolved Predicted and Experimental Propfan 

Surface Static Pressure Histories for SR7 Propfan - Cylindrical Post 
Aerodynamic Interaction Study (36.7% Axial Chord, 64.1% Radial 
Span, Pressure Side) 96 




Figure 4.36: Predicted Instantaneous Propfan Airfoil Midspan Total Pressure Con- 

tours for Cylinder /Propfan Interaction Study 
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also described in Section 4.2. Geometric and aerodynamic design parameters for the 
model propeller are given in Figure 4.4. The calculations presented here were based 
on a cruise condition with a flight Mach number of 0.72, and advance ratio of 3.14, 
and forward and rearward blade setting angles of 57.4 and 60.0 degrees, respectively. 
This time-dependent flow calculation is the unsteady equivalent of the steady flow 
mixing plane calculations described in Section 4.2. 

Time-dependent flow performance predictions for the CRP-X1 propfan were ob- 
tained by using the time-dependent inter-blade row interpolation concept described 
in Section 3.6. The grid therefore consists of two blocks (one for each blade row) as a 
result of the simple blade count ratio (1:1). It should be noted that this grid system 
corresponds to standard configuration #7 as described in the ADPAC-AOACR Com- 
puter Program Users Manual [9]. The grid block sizes were 65x49x33 and 57x49x33, 
respectively. An illustration of the mesh system is given in Figure 4.5. The time- 
dependent calculation was initiated from the steady flow predictions obtained using 
the mixing plane concept. In this case, because of the simple blade count ratio, 
the same grid system may be used for either calculation. Once initiated, the time 
dependent solution updates the data between the blade row mesh blocks based on a 
time-space resolved interpolation procedure to accurately maintain the interblade row 
aerodynamic interactions. The solution eventually becomes time-periodic, as shown 
in the residual history given in Figure 4.37. An illustration of the time-dependent 
aerodynamic interactions which occur between the counterrotating propfan airfoils is 
illustrated in the instantaneous midspan static pressure contour plots given in Fig- 
ure 4.38. This series of instantaneous plots illustrates the relative motion between 
the propfan blade rows and the resulting static pressure field interactions which result 
from this motion over one cycle of the interaction period at the blade midspan. The 
interaction is obviously most intense at the point where the adjacent propfan airfoils 
are lined up circumferentially. The interaction appears to have a more significant 
effect on the downstream blade row, which is not unexpected since the downstream 
blade row is subject to both potential and vortical disturbances while the upstream 
blade row is only subject to potential disturbances.. In the absence of time-dependent 
experimental data, these results can only be interpreted philosophically; however, the 
overall indications are that the time- dependent flowfield has at least been captured 
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Model Counterrotating Propfan 



Figure 4.37: ADPAC-AOACR Time- Accurate Residual History for Model Counter- 

Rotating Propfan Aerodynamic Interaction Analysis 
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in a reasonable manner, and has been done relatively efficiently by taking advantage 
of the capabilities of the ADPAC-AOACR code. 


4.8 GMA3007 Fan Section Distortion Study (Unsteady Flow) 

The final calculation to be presented is a prediction of the unsteady flow which 
results from a circumferential inlet distortion in a modern turbofan engine fan sec- 
tion. This calculation serves the dual purpose of analyzing a complex aerodynamic 
phenomena associated with turbofan engine operation, as well as displaying the cou- 
pled 2-D/3-D mesh block capabilities of the ADPAC-AOACR code. The geometry 
selected for this study is the fan section of the Allison GMA3007 turbofan engine. 
The Allison GMA3007 is a 5:1 bypass ratio turbofan engine which produces 7,000 
pounds of thrust. The fan section geometry considered consists of three blade rows 
including the fan rotor, bypass stator, and core stator. The flow is complicated by the 
presence of the core flow splitter, and the fact that there are multiple exits (bypass 
flow, core flow) for which proper boundary data must be specified. An outline of the 
GMA3007 fane section geometry is given in Figure 4.17. 

The simulation selected involves a 6 per revolution circumferential distortion. 
The distortion is represented by a stationary total pressure deficit in the incoming 
flow. Each area of distortion occupies a 30 degree arc circumferentially, and is located 
between 60% and 100% span radially. This particular distortion pattern was chosen 
because the resulting mesh block structure can be simplified considerably compared 
to other (1,2, or 3 per revolution) distortion patterns. The radial distribution of the 
distortion was chosen based on experience which suggests that most fan designs are 
“tip critical”, suggesting that distortion effects are most critical near the fan tip. The 
distortion pattern is illustrated graphically in Figure 4.39. 

This calculation required a rather unique mesh block arrangement as shown in 
Figure 4.40. The fan section consists of three blade rows (fan rotor, fan bypass 
stator, and the core stator). The axial separation between blade rows suggest that it 
is possible to solve for the rotor only in this calculation; however, previous experience 
suggests that the presence of the downstream geometry, and particularly the core 
flow splitter may have some influence on the time-averaged aerodynamic loading 
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Figure 4.38: ADPAC-AOACR Predicted Instantaneous Static Pressure Contours at 

Midspan (1 cycle of time-periodic solution) 
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Figure 4.39: Total Pressure Deficit Pattern for Allison GMA3007 Engine Fan Section 

Distortion Study 
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of the upstream rotor. Four stationary mesh blocks located upstream of the fan 
rotor are utilized to impose the desired circumferential and radial distortion pattern. 
The distortion is defined as a nonuniform inlet boundary condition on two of the 
stationary mesh blocks, resulting in the distortion pattern previously described. The 
fan rotor is represented by a mesh system discretizing four blade passages. This is the 
minimum portion of the rotor which must be represented under the present distortion 
pattern. The use of four blade passages permits an application of periodicity across 
the rotor blade representation which is consistent with the imposed 6 per revolution 
distortion pattern. It was assumed that the bypass stator and core stator did not 
require more than a single blade passage representation, and therefore the bypass 
stator is represented by a single 3-D blade passage mesh, and the core stator is 
represented through a 2-D mesh system with embedded blade element blockage and 
body forces to simulate the axisymmetric effects of the stator. This representation 
was chosen to demonstrate the combined 2-D /3-D solution coupling capabilities of 
the ADPAC-AOACR code. The single passage mesh blocks for the bypass stator 
and the 2-D representation of the core stator are numerically coupled with the mesh 
blocks for the fan rotor by using a circumferential averaging procedure. Due to time 
constraints and computational limitations, a relatively sparse mesh was utilized in this 
calculation, and the results should therefore be interpreted only as a demonstration 
of this capability, and not conclusively representative of the truie flow. 

The time-dependent calculation was initialized from a steady flow calculation 
utilizing the mixing plane concept between blade rows. The time-dependent calcu- 
lation was performed over approximately 4 cycles of the distortion pattern, beyond 
which point a time-periodic solution was observed. 

During the course of this calculation, the distortion pattern was observed to 
become less well defined with axial distance from the inlet. This phenomenon was a 
result of the relatively coarse mesh used in this demonstration which unrealistically 
smeared the distortion pattern due to numerical dissipation. As a result, the impact 
of the distortion on {he rotating fan was minimized by the resulting smoothing of 
the distortion as it approached the fan. This implies that for practical calculations 
involving inlet distortions, a significant number of grid points must be utilized just 
on the basis of distortion convection alone. 
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Figure 4.40: ADPAC-AOACR 10-Block Mesh System for Allison GMA3007 Engine 

Fan Section Inlet Distortion Study 
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An illustration of the effects of the inlet distortion on the fan rotor are displayed 
on the surface static pressure contour plot given in Figure 4.41. The effects of the 
distortion are displayed clearly in the repeating pattern near the tip of the fan rotor 
blades. Total pressure contours on a mesh plane downstream of the fan leading edge 
at approximately 30% axial chord illustrated in Figure 4.42. This figure illustrates 
the resulting total pressure distortion pattern after entering the rotor passage. Again, 
the distortion continues to become less noticable with axial distance as a result of 
artifical dissipation. The influence of secondary flow and aerodynamic loading on the 
airfoils transform the shape of the distortion, but the overall extent of the distortion 
appears to be relatively unchanged. Without experimental data to verify the results 
of the analysis, no concrete conclusions can be drawn from this analysis; however, the 
potential for the capability to analyze complex turbomachinery flows resulting from 
inlet flow distortion has been demonstrated. 
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Figure 4.41: ADPAC-AOACR Instantaneous Predicted Rotor Surface Static Pres- 

sure Contours for GMA3007 Fan Operating Under Inlet Distortion 
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Figure 4 42- ADPAC-AOACR Predicted Instantaneous Total Pressure Contours at 
‘ ° ' ' 30 % Axial Chord for the GMA3007 Fan Operating Under Inlet Distor- 

tion 
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5. CONCLUSIONS 


A time-dependent, three-dimensional Euler /Navier- Stokes aerodynamic analysis 
has been developed for the numerical analysis of ducted and unducted fan flowfields 
involving multiple blade rows. The underlying multi-block discretization scheme 
provides a convenient basis upon which detailed aerodynamic analyses for compli- 
cated counterrotating geometries may be performed. Multiple blade row calcula- 
tions were demonstrated utilizing both the steady flow mixing plane concept and the 
time-resolved interpolation procedure for rotor/stator interaction calculations. Aero- 
dynamic predictions generated from the analysis were verified through comparisons 
with both steady state and time-dependent experimental data. The capability of ac- 
curately simulating the time-dependent aerodynamics in an advanced turbofan engine 
fan section with inlet distortion has been demonstrated. 

Several comments are in order concerning the various numerical techniques ap- 
plied in this study. It is apparent that the algebraic turbulence model may not be 
well suited for the complex vortical flows which can occur for modern turbomachmery 
blade designs. Future efforts may benefit from more detailed turbulence models devel- 
oped for complex 3-D flows. The time-accurate implicit residual smoothing algorithm 
is also a known source of error for unsteady flow predictions. Although these errors 
are suspected to be local in nature, additional studies dedicated to understanding 
the impact and limitations of this algorithm may be warranted. The overall accuracy 
of the analysis can be swayed by additional factors, including the unknown deflected 
shape of the propfan blade, errors introduced through poor grid resolution, turbulence 
modeling, and artificial dissipation. In spite of the known algorithmic deficiencies, the 
analysis has successfully predicted several complex time-dependent flows with inter- 
acting geometries, and has demonstrated good agreement with available experimental 

data. 
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